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I.  INRT ODCUCTION 


The  postdoctoral  fellowship  grant  was  awarded  to  the  principal  investigator  (PI)  for  the  period  of 
April  1,  2003 — March  31,  2005.  The  purpose  of  this  investigation  is  to  introduce  a  framework  for 
including  model  parameter  uncertainties  into  prostate  Intensity  Modulation  Radiation  Therapy 
(IMRT)  dose  optimization  so  that  biological  model-based  objective  function  can  be  used  with 
improved  confidence  level.  The  specific  aims  of  the  proposal  are:  (1)  to  establish  a  mathematical 
formalism  to  incorporate  model  parameter  uncertainty  into  IMRT  optimization;  (2)  to  identify  the 
clinically  relevant  biological  model  parameter  variance  range;  and  (3)  to  study  the  prostate  cancer 
treatment  planning  including  the  model  uncertainty  information.  Under  the  generous  support  from 
the  U.S.  Army  Medical  Research  and  Materiel  Command  (AMRMC),  the  PI  has  contributed 
significantly  to  the  radiation  treatment  of  prostate  cancer.  Several  conference  abstracts  and  refereed 
papers  have  been  resulted  from  the  support.  The  fellowship  also  allowed  the  PI  to  obtain  research 
training  in  prostate  cancer  while  accomplishing  the  proposed  projects.  The  preliminary  data  and 
research  opportunity  gained  under  the  support  of  this  grant  has  enabled  the  PI  to  obtain  offers  of  an 
assistant  professor  in  the  Department  of  Radiation  Oncology  at  a  few  prestigious  universities. 

II.  RESEARCH  AND  ACCOMPLISHMENTS 

Adenocarcinoma  of  the  prostate  is  the  most  common  malignancy  in  men  in  the  western 
countries.  Options  for  active  management  of  organ-confined  prostate  cancer  include  radical 
prostatectomy  and  definitive  radiotherapy  with  either  external  beams  or  interstitial  brachytherapy. 
Intensity  Modulated  Radiation  Therapy  (IMRT)  is  quickly  replacing  conventional  techniques  for  the 
treatment  of  prostate  cancer.  Most  IMRT  optimization  systems  at  present  use  dose  and/or  dose 
volume-based  objective  functions1,  which  guide  the  IMRT  planning  by  imposing  a  penalty 
according  to  the  difference  between  the  computed  and  prescribed  doses.  A  well-known  drawback  of 
the  dose-based  inverse  planning  is  that  the  nonlinear  dose  response  of  tumor  or  normal  structures  is 
not  fully  considered.  A  number  of  mathematical  models  have  been  developed  over  the  years  to 
better  describe  the  biological  effect  of  radiation  and  considerable  works  have  also  been  done  to  use 
these  biological  models  to  construct  more  meaningful  objective  functions  for  therapeutic  dose 
optimization2.  Generally  speaking,  radiobiological  formalism  involves  the  use  of  model  parameters 
that  are  of  considerable  uncertainty3'6.  For  instance,  the  radiosensitivity  of  Webb’s  TCP  model 
varies  from  0.157  Gy'1  to  0.090  Gy"'  when  model  parameters  were  fit  to  103  patients’  data3.  In  order 
to  improve  the  dose  distribution  and  the  outcome  of  prostate  cancer  radiotherapy,  we  implemented  a 
mathematical  formalism  to  include  all  types  of  model  parameter  uncertainties  to  the  dose 
optimization  in  the  first  part  of  this  project.  Based  on  this  result,  in  second  and  third  parts  of  work, 
we  developed  a  biological  optimization  framework  and  incorporated  the  model  parameter 
uncertainties  into  IMRT  optimization  based  on  the  clinical  outcome  knowledge.  We  also 
implemented  an  algorithm  to  optimize  the  time-dose-fractionation  for  the  radiation  treatment  of 
prostate  cancer  with  inclusion  of  the  biological  model  parameter  uncertainties. 

An  important  issue  in  inverse  treatment  planning  is  how  to  formalize  the  clinical  goals  to 
objectively  evaluate  the  figures  of  merit  of  different  IMRT  plans.  Over  the  last  two  decades, 
attempts  have  been  made  by  many  researchers  to  capture  the  main  feature(s)  of  the  dose  volume 
effects.  A  power  law  model  represents  one  of  the  successful  techniques  in  dealing  with  the  dose- 
volume  effects  of  sensitive  structures7.  In  this  model  an  equivalent  dose  uniformly  irradiating  the 
whole  organ,  Deq,  can  be  used  to  represent  the  situation  in  which  a  fractional  partial  volume,  v,  is 
irradiated  to  a  dose,  D,  by  a  simple  power  law  model:  Deq=vI/nD.  A  remarkable  characteristic  of  this 
model  is  that,  although  only  a  single  organ-specific  parameter,  n,  is  used,  clinical  and  biological  data 
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has  shown  that  this  power  law  holds  well  at  low  complication  levels7, 8.  Based  on  this  relation, 
Mohan  et  al9  introduced  the  concept  of  effective  dose  to  represent  a  non-uniform  dose  distribution  in 
a  sensitive  structure.  Kutcher  and  Burman10  applied  the  same  power  model  independently  to  each 
volume  element  of  the  histogram  and  introduced  the  concept  of  effective  volume  to  reduce  the  DVH 
of  an  inhomogeneous  dose  distribution  in  a  sensitive  structure  to  a  uniform  dose  distribution. 
Following  their  study,  in  this  project  we  define  the  effective  volume  (A Vejf)i  for  a  voxel  i  with 
volume  AF and  dose  Dc(i)  as  follows 

(A  Veff)i  =  AV(DC  (/)  /  Dref)Un  (1) 


and  extend  this  concept  to  handle  the  voxels  in  the  tumor  target,  where  n  is  an  organ-dependent 
parameter  and  Dref  is  the  reference  dose.  For  a  sensitive  structure,  n  is  a  small  positive  number 
(0<n<l)  and  the  value  of  parameter  n  reflects  the  architecture  (serial  or  parallel)  of  the  sensitive 
structure.  For  a  target,  n  should  be  assigned  with  a  small  negative  value  (-1<«<0). 

The  objective  function,/,  expressed  as  a  function  of  the  effective  volume  in  the  voxel 
domain  for  an  organ  should  take  the  form  of 

/  =  /({(  AF*),}),  (2) 

A  more  general  form  of  inverse  planning  objective  function  can  be  written  as  a  hybrid  of  the 
dose-volume  based  and  the  dose-based  functions.  In  this  situation,  the  overall  objective  function  of 
the  system  takes  the  form  of 

f  -  i>,  /£<• + •)/ d,,„  ]■'•■  }|d,(o  -  di  (,f 

T=\  j=| 


+ 2X  /  £  o + <[^(0 1  d.m  V"")d<  (o'- . 


(7=1 


/=! 


(3) 


where  tt  and  sa  are  the  numbers  of  targets  and  sensitive  structures,  D'0  (/)  is  the  prescription  dose  in 


target  voxel  i,  subscripts  r  and  a  represent  target  r  and  sensitive  structure  a,  Nr,  Na,  rT,  ra  nT,  na,  Dr,ref, 
Da, ref,  K,  and  k„  represent  the  total  numbers  of  voxels,  structure  specific  importance  factors,  n 
parameters,  reference  doses,  power  of  dosimetric  deviation  from  the  specified  criteria  for  target  t 

and  sensitive 

structure  cr, 

respectively.  The 
factor 

Dc{i)-DT0(if  for 
target  or  Dc(i)k°  for 


(b)  Axial  1  Axial  2  Sagittal 

Fig.  1.  Comparison  of  the  isodose  distributions  of  the  two  prostate  IMRT  plans:  (a)  the 
conventional  dose-based  approach;  (b)  the  newly  proposed  approach.  The  results  on  two 
transverse  slices  and  a  sagittal  slice  are  shown. 


a  sensitive  structure 
represents  the 

contribution  from 
dosimetric  deviation 
from  the  ideal 
situation.  If  the  kx 
and  ka  are  set  to 
zero,  the  objective 
function  becomes 
purely  dose-volume 
driven.  In  particular, 
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n 


if  we  set  k„  to  zero  and  kr  to  a  non-zero  value,  the  objective  function  for  a  target  becomes  a  hybrid  of 
dose-volume  and  dose  based,  whereas  the  objective  functions  for  critical  structures  remain  to  be 
purely  dose-volume  based.  On  the  other  hand,  when  all  the  n  parameters  in  Eq.  (4)  are  set  to  be  +oo , 
no  dose-volume  effects  are  considered  and  Eq.  (4)  is  reduced  to  the  conventional  dose-based 
objective  function. 


Fig.  2.  Comparison  of  Dose  Volume 
Histograms  (DVHs)  of  the  prostate  IMRT  plans 
obtained  using  the  proposed  approach  (solid 
curves)  and  the  conventional  dose-based 
approach  (dash  lines). 


The  proposed  optimization  algorithm  is  used 
to  study  30  prostate  cancer  cases  and  the  results 
were  compared  with  that  of  the  conventional  dose 
based  IMRT  optimization  technique.  Figures  1  and  2 
are  the  results  of  the  two  IMRT  plans  obtained  using 
the  proposed  and  conventional  techniques  for  a 
typical  prostate  case.  Figure  1  compares  the  isodose 
distributions  in  two  transverse  slices  and  a  sagittal 
slice  for  the  two  plans.  The  DVHs  of  the  target  and 
sensitive  structures  are  plotted  in  figure  2,  in  which 
the  solid  and  dashed  lines  represent  the  DVHs 
obtained  using  the  new  and  conventional 
approaches,  respectively.  It  is  found  that,  for 
comparable  target  coverage,  the  new  inverse 
planning  technique  greatly  improves  the  critical 
structure  sparing,  especially  the  rectum  sparing. 
Furthermore,  it  is  intriguing  that  the  non-sensitive 
structure  normal  tissue  also  receives  fewer  doses  in 
comparison  with  that  of  the  dose-based  optimization. 
Our  results  suggest  that  the  improvement  in  the 
critical  structure  sparing  is  achieved  not  at  the  cost 
of  higher  target  dose  inhomogeneity,  which  is 
commonly  seen  in  IMRT  plan  optimization.  The 
calculated  NTCPs  of  rectum,  bladder  and  femoral 
heads  for  both  IMRT  plans  are  listed  in  table  I. 
According  to  the  table,  it  is  seen  that  the  NTCPs  of 
the  sensitive  structures  are  improved  significantly. 


For  the  rectum,  for  example,  the  Table  I  Comparison  of  the  normal  tissue  complication  probabilities 
NTCP  is  reduced  from  0.45%  to  (NTCP)  for  the  two  IMRT  plans  for  the  prostate  case 


0.03%. 

We  also 

implemented  a  method 
for  optimizing  the 
treatment  protocols  for 
prostate  cancer  with  the 

biological  model  parameter  uncertainties  included.  The  study  starts  from  an  extended  LQ  model 
with  inclusion  of  the  “4-R’s”11,  12 .  The  optimum  dose-time-fractionation  is  then  formulated  as  a 
problem  of  searching  for  the  highest  tumor  biologically  effective  dose  (BED)  while  keeping  the 
normal  tissue  BED  constant.  A  salient  feature  of  the  technique  is  that  various  influencing 
radiobiology  parameters,  such  as  the  redistribution  and  reoxygenation,  are  incorporated  naturally 


NTCP  (%) 

The  dose-based  IMRT  plan 

The  proposed  IMRT  plan 

Bladder 

0.017 

0.00030 

Rectum 

0.45 

0.029 

Femoral  head  (R) 

0.000076 

0.0000038 

Femoral  head  (L) 

0.000032 

0.000015 
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during  the  optimization  process.  The  optimized  tumor  BED  as  a  function  of  the  overall  treatment 
time  for  different  potential  doubling  time,  (15,  30  and  40  days),  is  shown  in  figure  3.  It  is  shown 


Figure  3.  The  optimized  tumor  BED  as  a  function  of  the 
overall  treatment  time  when  F  .  =15, 30,  40  days. 


that  the  overall  treatment  time  should  be 
larger  than  a  minimum  value  to  maximize 
the  tumor  BED  when  the  resensitization 
effect  is  considered  and  the  minimum 
overall  time  slightly  depends  on  the 
potential  double  time  Tpd- 

The  optimum  fractional  dose 
distributions  with  the  maximum  fractional 
dose  constraints  set  to  3  Gy  and  5  Gy  are 
shown  in  figure  4.  In  figures  4a  and  4b, 
the  a/|3  ratio  for  tumor  is  1 .5  Gy  and  the 
overall  time  is  43  days.  In  figures  4c  and 
4d,  the  a/p  ratio  for  tumor  is  3.0  Gy  and 
the  overall  time  is  38  days.  It  is 
interesting  to  observe  that  many  fractional 
doses  become  zero  and  a 
hypofractionation  scheme  with  the  size  of 


the  maximum  fractional 
dose  constraints  is  more 
favorable.  The  non-zero 
fractionations  are 
almost  equally  spaced 
over  the  entire 
treatment  time  and  the 
optimum  number  of 
fractionation  is 

determined  mainly  by 
the  maximum  fractional 
dose  constraint.  For 
example,  the  number  of 
optimum  fractionation 
is  20  and  9, 
respectively,  for  the 
maximum  fractional 
dose  constraint  of  3  and 
5  Gy.  Our  results 
indicate  that 

hypofractionation 
remains  to  be  the 
optimum  treatment 
scheme  even  when  the 
a/p  ratio  for  tumor  is 
the  same  as  that  of  the 
late  responding  tissue 
(both  are  3.0  Gy). 
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Figure  4.  The  optimum  fractionation  doses  with  ts  =  2  days  for  slowly  proliferating 
tumors,  (a)  a/p  =  1.5  Gy,  the  maximum  fractional  dose  constraint  is  3  Gy  and  the  overall 
treatment  time  is  43  days;  (b)  a/p  =  1.5  Gy,  the  maximum  fractional  dose  constraint  is  5 
Gy  and  the  overall  treatment  time  is  43  days;  (c)  a/p  =  3  Gy,  the  maximum  fractional 
dose  constraint  is  3  Gy  and  the  overall  treatment  time  is  38  days;  and  (d)  a/p  =  3  Gy,  the 
maximum  fractional  dose  constraint  is  5  Gy  and  the  overall  treatment  time  is  38  days. 
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III.  KEY  RESEARCH  ACCOMPLISHMENTS 


•  Developed  a  series  of  mathematical  formulae  to  incorporate  model  parameter  uncertainty 
into  IMRT  optimization. 

•  Verified  the  new  method  in  the  prostate  cancer  case  and  better  dose  coverage/sparing  can  be 
obtained  with  appropriate  parameters. 

•  Developed  a  clinical  knowledge-based  IMRT  inverse  treatment  planning  with  the  obtained 
model  parameters  and  studied  30  prostate  cancer  cases  using  the  proposed  inverse  planning 
framework. 

•  Implemented  an  algorithm  to  optimize  the  time-dose-fractionation  for  the  radiation  treatment 
of  prostate  cancer  with  inclusion  of  the  biological  model  parameter  uncertainties. 

IV.  REPORTABLE  OUTCOMES 

The  following  is  a  list  of  publications  resulted  from  the  grant  support.  Copies  of  the  publication 
materials  are  enclosed  with  this  report. 

Refereed  publication: 

1.  Yang  Y.  and  Xing  L.  Optimization  of  radiotherapy  dose-time-fractionation  with  consideration  of 
tumor  specific  biology.  International  Journal  of  Radiation  Oncology  Biology  and  Physics, 
(submitted),  2005. 

2.  Yang  Y  and  Xing  L.  Towards  Biologically  Conformal  Radiation  Therapy  (BCRT):  Selective 
IMRT  Dose  Escalation  Under  the  Guidance  of  Spatial  Biology  Distribution.  Medical  physics* 
(Accepted),  2005. 

3.  Yang  Y.  and  Xing  L.  Clinical  knowledge-based  inverse  treatment  planning.  Physics  in  Medicine 
and  Biology  49:  5101-51 17(2004) 

4.  Lian  J.  and  Xing  L.  Incorporating  Model  Parameter  Uncertainty  into  Inverse  Treatment  Planning, 
Medical  Physics  31:  271 1-2720  (2004). 

5.  Lian  J.,  Cotrutz  C.  and  Xing  L.  Therapeutic  treatment  plan  optimization  with  probability  density- 
based  dose  prescription.  Medical  Physics  30:  655-666  (2003). 

Published  Abstracts: 

1.  Yang  Y.  and  Xing  L.  Inverse  Treatment  Planning  with  Adaptively  Determined  Voxel-Dependent 
Importance  Factor,  46th  annua!  Meeting  of  AAPM,  Pittsburgh,  July  2004. 

2.  Yang  Y.  and  Xing  L.  Clinical  knowledge-based  inverse  treatment  planning.  In:  the  46th  American 
Society  for  Therapeutic  Radiation  and  Oncology  (ASTRO)  meeting,  Atlanta,  2004. 

3.  Tan  C.,  Yang  Y.,  Boyer  A.  and  Xing  L.  Enhancing  the  efficacy  of  radiation  therapy  by 
incorporating  spatial  distribution  of  heterogeneous  biology.  In:  the  46th  American  Society  for 
Therapeutic  Radiation  and  Oncology  (ASTRO)  meeting,  Atlanta,  2004. 

4.  Lian  J.,  Spielman  D.,  Cotrutz  C.,  Hunjan  S.,  Adalsteinsson  E.,  King  C.,  Luxton  G.,  Kim  D., 
Daniel  B.  and  Xing  L.  Including  metabolic  uncertainty  into  proton  MR  spectroscopic  imaging 
(MRSI)-guided  inverse  treatment  planning.  In:  the  45th  American  Association  of  Physicists  in 
Medicine  (AAPM)  Meeting,  San  Diego,  2003 

5.  Lian  J.  and  Xing  L.  Biological  Model  Based  IMRT  Optimization  with  Inclusion  of  Parameter 
Uncertainty.  In:  the  14th  International  Conference  on  the  Use  of  Computers  in  Radiation  Therapy 
(ICCR),  Seoul,  Korea,  2004 
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New  Employment: 

With  the  help  of  this  grant  and  associated  publications,  the  previous  PI  (Jun  Lian)  has  obtained  an 
assistant  professor  position  in  the  Department  of  Radiation  Oncology  at  the  University  of  North 
Carolina  at  Chapel  Hill.  The  Current  PI  (Yong  Yang)  has  obtained  offers  of  assistant  professor  in 
the  department  of  radiation  oncology  at  several  prestigious  universities  (written  offers  from  UCSF 
and  University  of  Pittsburgh,  and  verbal  offer  from  UCLA).  He  is  in  the  process  of  choosing  one  of 
them  to  start  his  academic  research  career.  He  plans  to  continue  his  research  in  prostate  cancer  to 
contribute  to  better  prostate  patient  care. 

V.  CONCLUSIONS 

Inverse  planning  is  an  important  step  in  IMRT  and  its  performance  crucially  determines  the 
quality  of  IMRT  treatment  plans.  In  this  work,  a  technique  for  incorporating  biological  model 
parameter  uncertainties  into  inverse  treatment  planning  has  been  developed.  By  including  model 
parameter  uncertainties,  we  provide  a  mechanism  for  incorporating  clinical  end  point  data  into 
inverse  treatment  planning  process  and  established  a  clinically  practicable  inverse  planning 
framework.  The  new  formalism  sheds  important  insight  into  the  problem  of  therapeutic  plan 
optimization.  The  results  of  30  prostate  cancer  cases  demonstrated  that  the  proposed  technique  is 
capable  of  greatly  improving  the  sensitive  structure  sparing  with  comparable  target  dose  coverage 
and  homogeneity.  In  addition,  through  accounting  for  the  known  uncertainties  in  the  model 
parameters,  we  implemented  an  algorithm  to  optimize  the  time-dose-fractionation  for  the  radiation 
treatment  of  prostate  cancer.  The  investigation  sheds  useful  insight  into  the  complex  dose-time- 
fractionation  problem  in  prostate  cancer  radiation  therapy  and  is  valuable  for  drafting  the  optimum 
clinical  trials  for  prostate  cancer  radiotherapy  and  for  interpreting  clinical  outcome  data. 
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Abstract 

Purpose:  To  explore  the  influence  of  the  “four  Rs”  of  radiobiology  on  external  beam 
radiotherapy  for  fast  and  slowly  proliferating  tumors  and  develop  an  optimization 
framework  for  tumor-biology  specific  dose-time-fractionation  scheme. 

Materials  and  Methods:  The  LQR  model  proposed  by  Brenner  et  al  (IJROBP,  32(2), 
1995)  is  used  to  describe  radiation  response  of  tumor,  in  which  the  time  dependence  of 
sublethal  damage  repair  is  included  and  redistribution  and  reoxygenation  effects  are 
described  using  a  term  of  resensitization  with  an  average  resensitization  time.  The 
optimum  radiotherapeutic  strategy  is  defined  as  the  treatment  scheme  that  maximizes 
tumor  biologically  effective  dose  (BED)  while  keeping  normal  tissue  BED  constant. 
Simulated  annealing  optimization  technique  is  used  to  search  for  the  optimal 
radiotherapeutic  strategies.  The  influence  of  different  model  parameters  on  total  dose, 
overall  treatment  time,  fraction  size  and  intervals  is  also  studied. 

Results:  For  fast  proliferating  tumors  the  optimum  overall  time  is  similar  to  the  assumed 
Tk,  the  time  from  the  beginning  of  treatment  to  the  starting  of  accelerated  proliferation, 
and  almost  independent  of  interval  patterns.  Significant  increase  in  tumor  control  can  be 
achieved  using  accelerated  schemes  for  the  tumors  with  doubling  time  smaller  than  3 
days,  but  little  is  gained  for  those  with  doubling  time  greater  than  5  days.  It  is  also  found 
that  the  incomplete  repair  of  normal  tissues  between  two  consecutive  fractions  in 
standard  fractionation  has  almost  no  influence  on  the  fractional  doses,  even  for  the 
hyperfractionation  with  an  interval  time  of  8h.  When  the  resensitization  effect  is 
included,  the  fractional  doses  at  the  beginning  and  end  of  each  irradiated  week  become 
obviously  higher  than  others  in  the  optimum  scheme  and  the  hyperfractionation  scheme 
has  little  advantage  over  the  standard  or  hypofractionation  one.  For  slowly  proliferating 
tumors,  provided  that  the  a/p  ratio  of  the  tumor  is  comparable  to  that  of  the  normal 
tissues,  a  hypofractionation  is  more  favorable.  The  overall  treatment  time  should  be 
larger  than  a  minimum,  which  is  predominantly  determined  by  the  resensitization  time. 
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Conclusion:  The  “four  Rs”  of  radiobiology  play  an  important  role  in  the  design  of 
radiation  therapy  treatment  protocol.  To  maximize  the  efficacy  of  radiation  therapy,  the 
properties  of  the  different  types  of  tissues  as  characterized  by  the  “four  Rs”  should  be 
exploited  and  incorporated  into  the  patient  treatment  through  the  optimization  of  dose- 
time-fractionation.  The  proposed  technique  provides  a  useful  tool  to  systematically 
optimize  radiotherapy  for  fast  and  slow  proliferating  tumors.  The  study  sheds  important 
insight  into  the  complex  problem  of  dose-time-fractionation  and  suggests  that  tumor  site- 
specific  optimization  has  great  potential  to  improve  therapeutic  outcome. 

Key  word:  Optimization,  linear-quadratic  model,  fractionation,  biological  model, 
radiobiology,  TCP,  NTCP. 
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INTRODUCTION 


Generally,  there  are  three  main  directions  to  improve  radiotherapy:  1)  improving  physical 
dose  distributions(l-6);  2)  optimizing  radiotherapy  dose-time-fractionation  scheme(7-10); 
and  3)  modifying  radiation  response  in  tumor  (radiosensitizers)(l  1,  12)  and/or  normal 
tissues  (radioprotectors)(13,  14)  using  chemical  agents.  With  clinical  implementation  of 
new  radiotherapy  techniques,  such  as  intensity-modulated  radiotherapy  (IMRT)(l-5),  it  is 
now  possible  to  achieve  optimal  physical  dose  distributions  on  a  patient  specific  basis. 
On  the  other  hand,  the  application  of  radiobiology  studies  in  clinical  practice  is  still  at  a 
primitive  stage  and  a  forward  comparison  of  tumor  control  probability  (TCP)  and  normal 
tissue  complication  probability  (NTCP)  is  often  employed  in  the  design  of  patient 
treatment  protocol.  To  facilitate  the  design  of  optimum  patient  treatment  protocols  and 
fully  exploit  the  potential  of  radiobiology  research  carried  out  over  the  years,  it  is  highly 
desirable  to  develop  a  technique  for  inverse  optimization  of  dose-time-fractionation  (15- 
20).  This  will  not  only  allow  us  to  better  understand  the  influence  of  various  biological 
parameters  on  radiation  therapy  treatment  of  cancer,  but  also  provide  a  practical  tool  to 
obtain  the  optimal  total  dose,  overall  treatment  time,  fraction  sizes  and  intervals  for  each 
disease  site. 

It  is  well  known  that  the  radiation  response  of  a  tissue  can  be  characterized  by  the 
“four  Rs”(21)  of  radiobiology:  repair  of  sublethal  damage,  repopulation,  redistribution 
and  reoxygenation.  Different  approaches  have  been  used  to  model  the  radiation  response, 
but  the  most  commonly  accepted  one  is  the  linear-quadratic  (LQ)  model(7,  22).  The 
original  LQ  model  describes  the  mechanism  of  cell  killing  and  captures  the 
characteristics  of  sublethal  damage  repair.  The  model  has  been  extended  to  include  time 
effect  and  the  influence  of  the  “four  Rs”(23,  24).  Current  standard  fractionation  and 
hyperfractionation  schemes  exploit  the  difference  in  the  repair  capability  to  sublethal 
damage  between  early-responding  tissue  and  late-responding  tissue  while  various 
accelerated  schemes  are  motivated  by  the  attempt  to  reduce  the  tumor  repopulation  effect. 
A  recent  research(9)  has  taken  the  difference  in  repair  rate  of  sublethal  damage  of  early- 
responding  and  late-responding  tissues  into  account  in  designing  brachytherapy 
protocols.  The  result  suggested  that  larger  doses  should  be  given  at  the  beginning  and  end 
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of  treatment  for  accelerated  treatment  regimens.  For  external  beam  radiotherapy, 
however,  little  systematic  study  has  been  done  to  exploit  the  temporal  processes  of  the 
repair,  repopulation  and  resensitization  (redistribution  and  reoxygenation)  and  their 
influence  on  the  optimum  treatment  strategy  for  different  types  of  tumors. 

In  this  work,  we  describe  a  method  for  optimizing  the  treatment  protocols  for  fast 
proliferating  tumors  (e.g.,  head  and  neck  cancer)  and  slowly  proliferating  tumors  (e.g., 
prostate  cancer).  The  study  starts  from  an  extended  LQ  model  with  inclusion  of  the  “4- 
R’s”.  The  optimum  dose-time-fractionation  is  then  formulated  as  a  problem  of  searching 
for  the  highest  tumor  biologically  effective  dose  (BED)  while  keeping  the  normal  tissue 
BED  constant.  A  salient  feature  of  the  technique  is  that  various  influencing  radiobiology 
parameters,  such  as  the  redistribution  and  reoxygenation,  are  incorporated  naturally 
during  the  optimization  process.  The  general  reference  drawn  from  this  study  is  that 
tumor  and  normal  tissue  biology  plays  a  significant  role  in  the  success  of  radiotherapy 
and  a  truly  individualized  treatment  is  highly  desirable  to  maximize  the  efficacy  of 
radiation  therapy. 


METHODS  AND  MATERIALS 

The  LQR  model  for  fractionation  scheme 

The  LQR  model  proposed  by  Brenner  et  al(24)  is  used  in  this  study  to  describe  radiation 
response  of  tumor.  In  this  model  the  time  dependence  of  sublethal  damage  repair  is 
included  and  the  redistribution  and  reoxygenation  effects  are  cast  in  a  term  of 
resensitization  with  an  average  resensitization  time.  The  surviving  fraction  S  of  tumor 
cells  irradiated  with  an  arbitrary  fractionation  scheme  can  be  expressed  as 

S  =  exp[  -aD-  fG(zR)D2  +(^ct2)G(ts)D2  +  H(T,al,Tk)(Tlo,  -Tk)ITfd],  (1) 

where  a  and  |3  are  linear  quadratic  constants  characterizing  the  intrinsic  radiosensitivity, 
<x2  is  variance  of  Gaussian  distribution  of  a,  is  repair  time,  is  average 
resensitization  time,  D  is  total  dose,  Ttot  is  overall  treatment  time,  Tpd  is  average  potential 
doubling  time,  and  Tk  is  “kick-off’  time  of  accelerated  proliferation,  representing  the  time 
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from  the  beginning  of  treatment  to  the  starting  of  accelerated  proliferation.  H(x,x0 ) 
denotes  the  Heaviside  function  and  is  defined  as 


H(x,x0) 


if  0<x<x0 
if  x>x0' 


G(t)  is  the  generalized  Lea-Catcheside  function  defined  by 


(2) 


G(t)  =  (~t)  duR(u )  J  dwR(w)  exp[-(w  -  w)!r\ ,  (3) 

where  R(t)  is  the  dose  rate  function.  For  an  arbitrary  fractionation  scheme  with  fractional 
dose  distribution  {dj...,  dt,...},  and  fractional  intervals  {Ati...,  At/,...}  ( d )  represents  the 
z'th  fractional  dose  and  Att  represents  the  time  interval  between  fractions  i  and  i+1). 
Assuming  that  the  delivery  time  for  each  fraction  is  ignorable,  we  can  rewrite  equation 
(1)  using  equation  (3)  as  follows 

5  -  exp {£[-«</,  -  f3G, (t, )df  +~*2G, (t, )df  ]  +  H(Tlol , Tk )(Tlnl  -Tk)/Tpd},  (4) 

i=]  £ 

where  7  is  total  fraction  number.  Here  Gfr)  is  calculated  for  each  fraction  with  following 
form: 

^(r)  =  l4Z^flexP(-A7,/r).  (5) 

“i  j= 1  k=j 

For  normal  tissues,  we  only  consider  the  cell  killing  and  incomplete  repair  of  sublethal 
damage,  i.e.,  only  the  first  two  terms  in  equation  (4)  are  included. 


Optimization  formulation 

The  goal  here  is  to  maximize  tumor  BED  while  keeping  normal  tissue  BED  constant. 
BED  is  calculated  using  the  cell  survival  S  given  by  equation  (4)  according  to  (17,  22): 

BED  =-— ln(5).  (6) 

a 

TCP  and  NTCP  are  calculated  using  (19,  25) 

TCP  or  NTCP-  exp(-A^)  =  exp[-770  exp(-aBED)] ,  (7) 

where  No  is  the  total  number  of  tumor  or  normal  tissue  cells.  To  proceed,  it  is  convenient 
to  use  the  BED  of  the  standard  scheme  (70Gy,  2Gy/f,  lf/d,  5f/w)  as  a  reference  because 
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its  properties  are  well  understood  and  extensive  clinical  data  for  the  scheme  has  been 
accumulated.  Mathematically,  the  dose-time-fractionation  optimization  can  be  found  by 
minimizing  the  objective  function 

F  -  BEDl0  /  BED, ,  (8) 

under  the  constraints 

BEDnl  =  BEDnl0  ,  (9) 

and 

BEDne<BEDne  o,  (10) 

where  BEDt,  BEDnU  and  BEDne  are  the  BEDs  of  the  target,  late-responding  normal  tissue 
and  early-responding  tissue  for  the  optimal  scheme,  respectively,  and  BEDl0,  BEDnl0 ?  and 
BEDne0  represent  the  same  for  the  standard  scheme. 

Theoretically,  the  interval  time  between  two  consecutive  fractions  can  be  treated 
as  one  of  the  system  variables  and  optimized  by  computer.  With  clinical  practicality  in 
mind,  we  only  study  two  patterns  of  intervals:  five  fractions  per  week  (interval=l  day) 
with  weekend  break  (standard  fractionation)  and  10  fractions  per  week  (interval  between 
two  fractions  in  the  same  day  is  8  h)  with  the  overnight  and  weekend  breaks 
(hyperfractionation).  The  overall  time  is  from  2  weeks  to  8  weeks  with  an  increment  of  1 
day.  For  each  overall  treatment  time,  our  optimization  program  will  then  search  for  the 
optimal  fractional  doses  that  maximize  the  tumor  BED  with  constant  normal  tissue  BED. 

By  introducing  Lagrange  multipliers,  Ai  and  Ae,  the  optimization  problem  with 
constraints(26,  27)  specified  in  equations  (9)  and  (10)  can  be  expressed  as 

F  =  BEDl0  /  BED,  +  A,  (. BEDnl  -  BEDnl0 ) 

+  AeH(BEDne ,  BEDne,  )(BEDne  -  BEDne0  ) ,  (11) 

where  H(x,x0 )  is  the  Heaviside  function  defined  in  equation  (2).  A  simulated  annealing 
optimization  technique(28,  29)  is  used  to  optimize  the  objective  function  (11)  with  the 
fractional  doses  as  the  optimization  variables.  For  each  given  interval  pattern,  the  optimal 
overall  treatment  time  is  found  by  comparing  the  optimal  tumor  BEDs  obtained  for  all 
possible  overall  treatment  times. 
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Estimation  of  Radiobiologic  Parameters 

Before  optimizing  the  dose-time-fractionation,  it  is  important  to  determine  the  values  or 
the  variation  ranges  of  the  model  parameters  for  the  studied  tumor  sites.  Seven 
independent  parameters  for  each  tumor  site  and  three  parameters  for  each  normal  tissue 


are  needed,  which  are  a,  a/p  ratio, 


Tk?  Tpd ,  ?R,  rs,  and  cr 


for  tumor  and  a,  a/p  ratio,  and 


zR  for  normal  tissue.  In  this  study  a  is  set  to  be  0.315/Gy  for  both  late-responding  and 
early-responding  normal  tissues,  a/p  ratio  and  repair  time  rR  are  chosen  to  be  3  Gy,  4  h 

and  10  Gy,  0.5  h  for  the  late-responding  and  early-responding  tissues,  respectively(10, 
22,  30).  For  fast  proliferating  tumors,  the  estimated  doubling  time,  Tpcp  is  around  3  days 
as  indicated  by  recent  studies(31,  32)  and  the  commonly  accepted  a  is  0.35/Gy  with  a/p 
ratio  =  10  Gy(10,  22).  Tk  for  head  and  neck  (FIN)  tumors  as  revealed  by  some  studies  is  in 

the  range  of  21  to  35  days(33),  and  the  repair  time  of  the  tumor  cell  rR  is  around  0.5  h(9, 
30).  The  resensitization  process  includes  both  the  redistribution  and  reoxygenation 
effects.  It  is  estimated  that  the  resensitization  time  ts  is  comparable  to  the  cell  cycle  time 


or  reoxygenation  characteristic  time,  i.e.,  hours  to  days.  In  this  work  we  assume  ts  to  be 


1  , 

in  the  range  of  0.5-3  days  for  fast  proliferating  tumors(24).  —a  represents  the  amplitude 


for  the  resensitization  and  is  assumed  to  be  0.02.  In  addition,  in  order  to  convert  the 
resultant  BED  to  TCP  we  assume  that  the  total  number  of  tumor  cells  No  =108. 

For  slowly  proliferating  tumors,  such  as  prostate  tumor,  a  is  assumed  to  be 
0.10/Gy  and  a/p  ratio  in  the  range  of  1. 0-4.0  Gy,  as  suggested  by  several  recent 
studies(19,  34,  35).  The  estimated  doubling  time  during  radiotherapy  is  around  40  days, 
with  a  range  from  16  to  61  days(36).  Because  there  is  no  Tk  data  available  for  slowly 
proliferating  tumors,  in  this  work  we  study  the  extreme  situation,  i.e.,  Tk  =0.  The  results 
should  be  extendable  to  situations  with  nonzero  Tk  values.  The  repair  time  of  prostate 
tumor  cells  is  1.9  h,  suggested  by  Fowler  et  al(35).  A  recent  study(37)  showed  that  there 
might  be  no  full  redistribution  of  cells  around  the  cell  cycle  between  radiation  fractions 

1  2 

spaced  by  1  day.  Thus,  we  assume  rs  to  be  in  the  range  of  1-3  days.  We  choose  —  cr  = 
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— (3  for  slowly  proliferating  tumors.  Table  1  summarizes  the  values  of  model  parameters 
used  in  this  work  for  the  fast  and  slow  proliferating  tumors. 


RESULTS 

Fast  proliferating  tumors 

Overall  treatment  time 

Figures  la  and  lb  show  the  optimized  tumor  BED  as  a  function  of  the  overall  treatment 
time  for  three  different  7)7  s  (21,  28  and  35  days)  for  the  two  studied  interval  patterns 
(standard  fractionation  and  hyperfractionation).  Other  parameters  used  to  obtain  the  data 
are  listed  in  Table  1.  It  is  found  that  there  exists  an  optimum  overall  treatment  time  with 
which  the  maximum  tumor  BED  is  achieved  for  each  Tk.  The  optimum  overall  time  is 
similar  to  the  assumed  Tk  and  almost  independent  of  interval  patterns.  For  the  standard 
fractionation,  for  example,  TCP  increases  by  about  60%  for  a  treatment  with  the  optimum 
overall  time  as  compared  with  the  conventional  7-week  overall  time  if  Tk  =28  days. 

Figures  2a  and  2b  show  the  BED  as  a  function  of  the  overall  treatment  time  for  a 
few  different  Tpd  values  (2,  3,  4  and  5  days)  for  the  two  studied  interval  patterns.  In  both 
cases,  the  influence  of  Tpd  is  more  pronounced  when  is  small.  When  TP*  is  larger  than 
5  days,  shortening  the  overall  time  results  in  almost  no  increase  in  tumor  BED.  For 
example,  the  difference  in  TCPs  between  the  treatment  schemes  with  the  optimum 
overall  time  and  the  conventional  7-week  overall  time  is  only  about  3%  when  days 
for  the  standard  fractionation. 

In  figures  3a  and  3b  we  show  the  results  for  Ts  =0.5,  1,  2,  to  3  days  for  the  two 
interval  patterns.  While  the  BED  is  sensitive  to  the  variation  of  rs,  the  optimum  overall 
time  is  not. 

Total  dose 
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The  optimized  total  dose  versus  the  over  time  for  a  few  Tk s  is  shown  in  figure  4  for  the 
two  interval  patterns.  As  expected,  the  total  dose  increases  with  the  overall  time. 
Generally,  the  fractional  dose  increases  with  the  reduction  of  overall  time.  In  order  to 
keep  normal  tissue  complications  (mainly  late  complications)  constant,  the  total  dose 
must  be  decreased.  Figure  4  also  suggests  that  the  hyperfractionation  require  a  higher 
total  dose  than  the  standard  fractionation  for  the  same  overall  time.  In  addition,  it  is  seen 
that  the  optimized  total  dose  is  independent  of  Tk  for  both  types  of  fractionation  regimes. 

Fractional  doses 

Figures  5  and  6  show  the  optimum  fractional  doses  obtained  for  two  examples:  one  for 
the  standard  fractionation  with  overall  time  of  40  days  and  total  fraction  number  of  29 
and  the  other  for  the  hyperfractionation  with  overall  time  of  46  days  and  total  fraction 
number  of  68.  In  figures  5a  and  6a  we  plot  the  results  with  xR  =0  h  for  both  tumor  and 
normal  tissues  and  rs  =0  h  for  tumor.  These  parameter  values  represent  the  situation  that 
the  sublethal  damage  repair  for  both  tumor  and  normal  tissues  and  the  resensitization  of 
tumor  cells  are  completed  immediately  after  irradiation.  In  this  situation,  it  is  found  that  a 
uniform  fractional  dose  is  more  favorable.  The  result  is  intuitively  conceivable  because 
evenly  distributed  fractional  doses  maximize  the  repair  of  the  late-responding  tissue, 
leading  to  the  maximum  tumor  BED  with  a  given  normal  tissue  damage.  The  influence  of 
different  normal  tissue  repair  time  on  the  optimum  fractional  doses  with  fixed  Ts  =0  h  and 

xR  =0.5  h  for  tumor  is  shown  in  figures  5b-c  and  6b-c.  For  the  standard  fractionation, 
when  late  tissue  repair  time  TR  is  set  to  be  4  h,  the  optimum  fractional  doses  remain 
uniform.  This  is  because  the  sublethal  damage  repair  of  normal  tissue  is  almost 
completed  during  the  interval  between  two  consecutive  fractions  (24  h).  When  tr  of 

normal  tissue  increases  from  4  h  to  12  h,  higher  fractional  doses  at  the  beginning  and  end 
of  each  week  become  preferable.  The  ratio  of  the  averaged  “spike”  dose  with  the 
averaged  dose  of  remaining  fractions  is  1.4.  For  the  hyperfractionation,  the  dose  “spikes” 
occur  even  if  the  normal  tissue  repair  time  rR  is  4  h.  The  ratio  in  this  case  is  1.1.  The 

effect  becomes  more  pronounced  as  xR  is  increased  to  8  h  and  the  ratio  becomes  1.4. 


10 


Figures  5d  and  6d  illustrate  the  results  for  the  average  tumor  resensitization  time  Ts  =24  h 
with  tR  =  0.5  h  for  tumor  and  4  h  for  the  late-responding  tissue.  Similar  to  the  situation 
discussed  above,  the  optimal  fractional  doses  at  the  beginning  and  end  of  each  week  are 
higher  when  the  resensitization  effect  is  included.  The  ratios  of  the  averaged  “spike”  dose 
with  the  averaged  dose  of  remaining  fractions  are  1.2  and  1.3  for  the  standard 
fractionation  and  hyperfractionation,  respectively. 


Interval  patterns 

A  comparison  of  the  tumor  BED  of  the  standard  fractionation  and  hyperfractionation  for 
the  same  late  tissue  complication  as  a  function  of  overall  treatment  time  is  plotted  in 
figure  7.  The  influences  of  different  normal  tissue  repair  time  ( xR )  and  tumor  average 


resensitization  time  (ts)  are  also  illustrated.  In  figure  7a  we  ignore  the  resensitization 
effect  and  assume  that  the  sublethal  damage  repair  is  completed  immediately  after 

irradiation  (i.e.,  -^-er2 ,  rs  and  rR  for  tumor  and  normal  tissues  are  all  set  to  be  zero).  The 


results  show  that,  in  this  situation,  the  BED  for  the  hyperfractionation  is  higher  than  that 
for  the  standard  fractionation.  Figures  7b  and  7c  show  the  results  with  the  incomplete 
repair  considered  but  the  resensitization  effect  ignored.  The  tumor  repair  time  is  set  to  0.5 
h  and  the  late  tissue  repair  time  4  h  and  8  h  in  figures  7b  and  7c,  respectively.  It  is  found 
that,  although  the  hyperfractionation  is  better  than  the  standard  fractionation,  the  BED 
difference  between  the  two  fractionation  regimes  decrease  as  the  late  tissue  repair  time 
increases.  Figure  7d  shows  the  results  with  the  resensitization  effect  included.  The 
average  resensitization  time  for  tumor  is  24  h  with  tr  =  0.5  h  for  tumor  and  4  h  for  the 
late  tissue.  As  seen  from  the  data,  the  BEDs  for  the  two  regimes  are  very  close  in  this 
situation. 


Slowly  proliferating  tumors 

When  the  a/p  ratio  for  tumor  is  lower  than  that  for  the  surrounding  normal  tissues,  our 
formulae  suggest  that  a  single  fractionation  irradiation  is  the  optimum  scheme  if  the 
resensitization  effect  is  ignored.  Clinically,  however,  a  single  fractionation  irradiation 
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may  result  in  disastrous  acute  sequelae  and  consequential  late  effect.  Because  there  is 
currently  no  good  way  to  model  the  acute  sequelae  and  consequential  late  effects,  in  this 
work  we  simply  restrict  the  fractional  dose  to  within  a  given  value  to  consider  these 
effects. 

Overall  treatment  time 

The  optimized  tumor  BED  as  a  function  of  the  overall  treatment  time  for  different 
(15,  30  and  40  days)  is  shown  in  figure  8a.  The  maximum  fractional  dose  is  constrained 
to  5  Gy  in  this  calculation.  Other  parameters  used  to  obtain  these  data  are  listed  in  Table 
1.  Different  from  fast  proliferating  tumors,  we  found  that  the  overall  treatment  time 
should  be  larger  than  a  minimum  value  to  maximize  the  tumor  BED  when  the 
resensitization  effect  is  considered.  The  optimum  BED  increases  gradually  with  the 
overall  time  before  this  minimum  overall  time  and  after  it  the  BED  almost  keeps 
constant.  It  is  also  shown  that  the  minimum  overall  time  slightly  depends  on  the  potential 
double  time  V  The  larger  the  Tpd'  the  larger  the  minimum  overall  time  is.  For  example, 
the  minimum  overall  time  is  about  5-6  weeks  for  Tp*  =  40  days,  but  about  3-4  weeks 
when  Tpd  =  1 5  days.  In  addition,  the  curve  for  Tpd  =  1 5  days  shows  that  the  tumor  BED 
slightly  decreases  after  about  6  weeks  due  to  the  influence  of  tumor  cell  proliferation. 
However,  the  decrease  disappears  if  Tk  takes  28  days  rather  than  0(data  not  shown). 

Figure  8b  gives  the  results  for  different  rs  (1,  2  and  3  days)  with  Tpd  =  40  days.  It  is  seen 
that  the  minimum  overall  time  is  significantly  dependent  on  rs:  the  larger  the  Ts,  the 
larger  the  minimum  overall  time.  The  minimum  overall  time  is  about  3  weeks  for  Ts  =  1 
day,  6  weeks  for  ts  =  2  days  and  larger  than  9  weeks  for  Ts  =  3  days.  This  conclusion 
implies  that  longer  overall  treatment  time  is  preferred  for  the  slow  proliferating  tumors 
(e.g,  tumors  with  more  hypoxic  cells). 

The  data  in  the  above  are  obtained  with  maximum  fractional  dose  constraint  of  5 
Gy.  The  calculation  can  be  easily  extended  to  other  constraint  values  and  similar 
conclusions  can  be  reached. 

Total  dose 
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Figure  9a  shows  the  optimized  total  dose  as  a  function  of  the  overall  treatment  time  when 
the  maximum  fractional  dose  is  3,  4,  5  and  6  Gy,  respectively.  The  parameters  used  to 
obtain  these  curves  are  listed  in  Table  1.  Unlike  the  case  of  fast  proliferating  tumors,  in 
this  circumstance,  total  dose  does  not  change  with  the  overall  treatment  time  because  the 
optimum  total  fractionation  number  is  purely  determined  by  the  given  maximum 
fractional  dose  constraint.  Furthermore,  as  shown  in  figure  9b,  total  dose  increases  with 
the  decrease  of  maximum  fractional  dose  constraints. 

Fractionation  scheme 

The  optimum  fractional  dose  distributions  with  the  maximum  fractional  dose  constraints 
set  to  3  Gy  and  5  Gy  are  shown  in  figure  10.  In  figures  10a  and  10b,  the  a/p  ratio  for 
tumor  is  1.5  Gy  and  the  overall  time  is  43  days.  In  figures  10c  and  lOd,  the  a/p  ratio  for 
tumor  is  3.0  Gy  and  the  overall  time  is  38  days.  It  is  interesting  to  observe  that  many 
fractional  doses  become  zero  and  a  hypofractionation  scheme  with  the  size  of  the 
maximum  fractional  dose  constraints  is  more  favorable.  The  non-zero  fractionations  are 
almost  equally  spaced  over  the  entire  treatment  time  and  the  optimum  number  of 
fractionation  is  determined  mainly  by  the  maximum  fractional  dose  constraint.  For 
example,  the  number  of  optimum  fractionation  is  20  and  9,  respectively,  for  the 
maximum  fractional  dose  constraint  of  3  and  5  Gy.  Our  results  also  indicate  that 
hypofractionation  remains  to  be  the  optimum  treatment  scheme  even  when  the  a/p  ratio 
for  tumor  is  the  same  as  that  of  the  late  responding  tissue  (both  are  3.0  Gy). 

The  change  of  the  optimized  tumor  BED  with  the  total  dose  with  tumor  a/p=1.5 
and  3  Gy  and  Ts=l  days  are  shown  in  figure  11a.  The  overall  treatment  time  is  46  days  for 
all  calculations  shown  in  figure  11a.  It  is  found  that  the  tumor  BED  decreases  with  the 
increase  of  total  dose  and  a  hypofractionation  schemes  yields  a  higher  tumor  BED  even  if 
the  tumor  a/p  ratio  is  the  same  as  that  of  the  late  responding  tissue.  Figure  1  lb  shows  the 
influence  of  different  ts  (1,  2,  and  3  days)  on  the  BED.  For  the  same  maximum  fractional 
dose  constraint,  the  optimized  tumor  BED  decreases  with  the  increase  of  resensitization 
time  Ts.  The  smaller  the  fractional  dose,  the  larger  the  influence  of  resensitization  time  Ts 
is.  Unlike  the  fast  proliferating  tumors,  the  fractional  dose  and  interval  time  do  not 
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change  with  Ts  as  shown  in  figures  1  lc  and  d,  in  which  the  Ts  is  set  to  1  day  and  3  days, 
respectively. 


DISCUSSIONS 

Several  clinical  and  experimental  studies  suggest  that  the  estimated  doubling  time  for  HN 
tumors  is  as  short  as  2-3  days(31,  32)  and  the  accelerated  tumor  growth  is  similar  for  all 
HN  tumor  sites  and  stages.  In  this  situation,  our  calculation  indicates  that  there  exists  an 
optimum  overall  treatment  time  that  maximizes  tumor  control  while  keeping  the  normal 
tissue  complication  constant.  We  found  that  the  optimum  overall  time  is  close  to  the 
assumed  Tk  and  almost  independent  of  fractionation  patterns  and  tumor  resensitization 
time,  which  is  similar  to  that  proposed  by  Fowler(lO).  However,  when  the  doubling  time 
is  larger  than  5  days,  it  is  found  that  the  tumor  control  improves  little  by  shortening  the 
overall  treatment  time,  suggesting  that  the  gain  of  accelerated  schemes  is  extremely 
sensitive  to  the  tumor  proliferation  rate.  The  conclusion  appears  to  be  consistent  with  the 
results  of  randomized  prospectively  controlled  clinical  trials  for  HN  tumors  by  the 

European  Cooperative  Radiotherapy  Group  (ECRG)  (38).  When  compared  with  the 

standard  fractionation  scheme,  the  ECRG  data  show  that  the  accelerated  scheme 

substantially  increases  the  local  control  rate  for  tumors  with  Tpd<  4  days,  whereas  for 

tumors  with  Tpd>  4  days,  there  is  no  detectable  difference  between  the  two  fractionation 
schemes.  Therefore,  in  order  to  truly  benefit  from  accelerated  fractionation,  it  is  essential 
to  select  patients  with  fast-proliferating  tumors. 

The  repair  time  is  typically  less  than  1  h  for  tumor  and  around  4  h  for  slow  repair 
component  of  normal  tissues,  as  suggested  by  a  number  of  clinical  and  laboratory-based 
animal  studies(9,  30,  39).  As  shown  in  this  work,  this  difference  in  the  repair  rate  has 
almost  no  influence  on  the  fractional  doses  for  the  standard  fractionation.  For 
hyperfractionation  with  two  consecutive  fractions  spaced  by  8  h,  the  difference  results  in 
slightly  higher  fractional  doses  at  the  beginning  and  end  of  each  week.  This  result  is  in 
accordance  with  the  dose  “spikes”  effect  mentioned  by  Brenner  et  al(9)  for  brachytherapy 
treatment  with  an  overall  treatment  time  of  120  h.  The  dose  “spikes”  effect  becomes 
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larger  with  the  lengthening  of  repair  time  for  normal  tissues  as  illustrated  in  this  study. 
The  underlying  biology  of  this  effect  is  that  normal  tissue  cells  can  get  better  repaired 
during  the  longer  rest  of  weekends.  Interestingly,  including  of  the  resensitization  effect 
leads  to  similar  results  as  shown  in  figures  5d  and  6d.  We  attribute  this  feature  to  the 
increased  tumor  cell  killing  caused  by  the  sufficient  resensitization  during  the  weekends. 
Furthermore,  this  work  indicates  that  the  studied  hyperfractionation  scheme  leads  to 
better  tumor  control  as  compared  to  the  standard  fractionation  scheme  even  when  the 
incomplete  repair  is  included.  However,  when  the  resensitization  effect  (ts=1  day)  is 

considered,  the  advantage  of  the  hyperfractionation  becomes  less  obvious  as  shown  in 
figure  7d,  suggesting  that  the  standard  fractionation  or  a  hypofractionation  scheme  may 
be  more  suitable  for  hypoxic  tumors. 

For  slowly  proliferating  tumors,  the  typical  doubling  time  is  about  40  days  and  the 
a/p  ratio  is  in  the  range  of  1  to  4  Gy.  Animal  studies  and  clinical  observations  support 
that  the  a/p  value  for  late  rectal  damage  is  greater  than  4  Gy(40).  Under  this 
circumstance,  it  is  obvious  that  hypofractionation  is  the  optimum  choice.  Our  work  also 
suggests  that  a  hypofractionation  scheme  is  preferable  even  if  a/p  ratios  for  both  tumors 
and  late  rectal  tissues  are  close,  e.g.,  all  equal  to  3.0  Gy. 

Different  from  fast  proliferating  tumors,  there  is  no  biological  advantage  to 
shorten  the  overall  treatment  time  for  tumors  with  slow  proliferation  rate.  Furthermore,  if 
the  resensitization  effect  is  considered,  there  exists  a  minimum  overall  treatment  time 
below  which  the  tumor  control  starts  to  decrease.  The  main  reason  responsible  for  this 
phenomenon  is  that  enough  interval  time  between  two  consecutive  fractions  is  required  to 
sufficiently  sensitize  the  tumor  cells.  Our  study  suggests  that  current  overall  treatment 
time  (6-7  weeks)  should  be  kept  for  the  hypofractionation  scheme.  This  will  enhance  the 
tumor  control  and  reduce  the  risk  of  acute  sequelae  and  consequential  late  effect  that  is 
not  easily  modeled  by  currently  available  radiobiology  theory.  In  addition,  as  shown  in 
figures  1  la  and  b,  although  a  larger  fractional  dose  results  in  a  higher  tumor  control,  the 
enhancement  becomes  insignificant  when  the  fractional  dose  is  greater  than  5-6  Gy  when 
a/p  ratios  for  the  tumor  and  later  rectal  tissue  are  comparable  or  when  the  resensitization 
time  is  large  (e.g.,  zs=3  days).  Therefore,  a  fractional  dose  of  5-6  Gy  seems  to  be  suitable 
for  prostate  hypofractionation  treatment.  We  note  that  the  risk  involved  in  this  kind  of 
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treatment  is  relatively  low.  The  experience  of  22  years  in  London  with  232  prostate 
cancer  patients  treated  with  a  6*6  Gy  scheme  in  3  weeks  has  indicated  comparable  local 
response  and  minimal  late  morbidity  when  compared  with  the  standard  fractionation 
scheme(41). 

Finally,  it  should  be  recognized  that  the  LQR  model  used  in  this  investigation  is 
far  from  satisfactory  in  modeling  the  responses  of  tumors  and  normal  tissues  and  there  is 
also  considerable  uncertainty  in  the  currently  available  model  parameters.  Notably,  using 
a  single  resensitization  time  for  both  redistribution  and  reoxygenation  effects  in  LQR 
model  should  be  improved  in  the  future  so  that  the  tumor  response  can  be  more  properly 
described.  In  order  to  yield  clinically  meaningful  results,  more  sophisticated  radiobiology 
model  and  more  accurate  parameters  are  needed.  In  addition,  theory  that  better  describes 
the  acute  reactions  and  consequential  late  effect  should  also  be  in  place. 


CONCLUSIONS 

In  this  work  we  have  described  a  general  framework  for  dose-time-fractionation 
optimization  and  explored  the  influences  of  the  “four  Rs”  of  radiobiology  on  radiation 
therapy  for  fast  and  slowly  proliferating  tumors.  Different  dose-time-fractionation 
schemes  are  evaluated  in  reference  to  the  standard  fractionation  and  the  role  of  various 
biological  parameters  in  the  design  of  a  treatment  protocol  is  elaborated.  Our  study 
indicates  that  it  is  clinically  significant  to  design  radiation  therapy  treatment  based  on  the 
specific  biological  properties  of  the  tumor  and  normal  tissues.  The  investigation  sheds 
useful  insight  into  the  complex  dose-time-fractionation  problem  in  radiation  therapy  and 
is  valuable  for  drafting  the  optimum  clinical  trials  for  different  tumor  sites  and  for 
interpreting  clinical  outcome  data. 
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Table  1  Model  parameters  used  in  our  study 


Parameters 

Fast  Proliferating 
Tumors 

Slowly  Proliferating 
Tumors 

LQ  constant  for  tumor,  a 

0.35/Gy 

0.10/Gy 

a/p  ratio  for  tumor 

lOGy 

1.5Gy 

“Kick-off’  time  for  tumor,  Tk 

28  days 

0  days 

Doubling  time  for  tumor,  Tpd 

3  days 

40  days 

Resensitization  time  for  tumor,  Ts 

1  days 

2  days 

Repair  time  for  tumor,  xR  tumor 

0.5  h 

1.9  h 

Variance  of  Gaussian  distribution 

-p 

,  1  2 

0.02 

of  a,  —a 

2 

3 

LQ  constant  for  late-responding 
normal  tissue,  a 

0.315/Gy 

0.315/Gy 

a/p  ratio  for  late-responding 
normal  tissue 

3  Gy 

3  Gy 

Repair  time  for  late-responding 

4  h 

4  h 

normal  tissue,  rRL  tissur 

LQ  constant  for  early-responding 
normal  tissue,  a 

0.315/Gy 

0.315/Gy 

a/p  ratio  for  early-responding 
normal  tissue 

10  Gy 

10  Gy 

Repair  time  for  early-responding 

0.5  h 

0.5  h 

normal  tissue,  rR  E_tissue 
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Legends 


Figure  1.  The  optimized  tumor  BED  as  a  function  of  overall  treatment  time  when  Tk  =21, 
28  and  35  days,  (a)  Results  for  standard  fractionation  (lf/d,  5  f/w);  and  (b)  Results  for 
hyperfractionation  (2f/d,  10  f/w). 

Figure.  2.  Influence  of  doubling  time  Tpd  (2,  3,  4  and  5  days)  on  the  optimum  overall 
treatment  time,  (a)  Results  for  standard  fractionation  (lf/d,  5  f/w);  and  (b)  Results  for 
hyperfractionation  (2f/d,  10  f/w). 

Figure  3.  Influences  of  average  resensitization  time  rs  (0.5,  1,  2,  and  3  days)  on  the 
optimum  overall  treatment  time.  Parameters  used  to  obtain  this  curves  are  listed  in  table 
1.  (a)  Results  for  standard  fractionation  (lf/d,  5  f/w);  (b)  Results  for  a  hyperfractionation 
pattern  (2f/d,  10  f/w). 

Figure  4.  The  optimized  total  dose  as  a  function  of  overall  treatment  time  when  Tk  =  21, 
28  and  35  days. 

Figure  5.  The  optimum  fractional  doses  for  four  different  sets  of  tR  and  Ts  for  the  standard 
fractionation.  The  overall  treatment  time  is  40  days  and  the  total  fraction  number  is  29. 

Figure  6.  The  optimum  fractionation  dose  for  four  different  sets  of  rR  and  rs  for  the 
hyperfractionation.  The  overall  treatment  time  is  46  days  and  total  fraction  number  is  68. 

Figure  7.  A  comparison  of  the  tumor  BEDs  under  the  same  normal  tissue  BED 
constraints  for  the  standard  fractionation  and  hyperfractionation  schemes. 

Figure  8.  The  optimized  tumor  BED  as  a  function  of  the  overall  treatment  time  when 
=15,  30,  40  days  (a)  and  ts=  1,  2,  3  days  (b).  The  maximum  fractional  dose  is  restricted 
to  5  Gy. 
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Figure  9.  (a)  The  optimized  total  doses  as  a  function  of  overall  treatment  time  for 
different  maximum  fractional  dose  constraints  (3,  4,  5  and  6  Gy),  (b)  The  optimized  total 
doses  as  a  function  of  the  overall  fractionation  number. 

Figure  10.  The  optimum  fractionation  doses  with  rs  =  2  days  for  slowly  proliferating 
tumors,  (a)  a/p  =  1.5  Gy,  the  maximum  fractional  dose  constraint  is  3  Gy  and  the  overall 
treatment  time  is  43  days;  (b)  a/p  =  1 .5  Gy,  the  maximum  fractional  dose  constraint  is  5 
Gy  and  the  overall  treatment  time  is  43  days;  (c)  a/p  =  3  Gy,  the  maximum  fractional 
dose  constraint  is  3  Gy  and  the  overall  treatment  time  is  38  days;  and  (d)  a/p  =  3  Gy,  the 
maximum  fractional  dose  constraint  is  5  Gy  and  the  overall  treatment  time  is  38  days. 

Figure  11.  (a)  The  optimized  tumor  BED  as  a  function  of  total  dose  for  different  a/p 
ratios  and  Ts.  (b)  The  optimized  tumor  BED  as  a  function  of  total  doses  for  different  rs; 
(c)  and  (d)  The  optimum  fractional  doses  for  rs  =  1  and  3  days.  In  all  calculation  the 
overall  time  is  set  as  46  days. 
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Abstract 


It  is  well  known  that  the  spatial  biology  distribution  (e.g.,  clonogen  density, 
radiosensitivity,  tumor  proliferation  rate,  functional  importance)  in  most  tumors  and 
sensitive  structures  is  heterogeneous.  Recent  progress  in  biological  imaging  is  making  the 
mapping  of  this  distribution  increasingly  possible.  The  purpose  of  this  work  is  to 
establish  a  theoretical  framework  to  quantitatively  incorporate  the  spatial  biology  data 
into  IMRT  inverse  planning.  In  order  to  implement  this,  we  first  derive  a  general  formula 
for  determining  the  desired  dose  to  each  tumor  voxel  for  a  known  biology  distribution  of 
the  tumor  based  on  a  linear-quadratic  (LQ)  model.  The  desired  target  dose  distribution  is 
then  used  as  the  prescription  for  inverse  planning.  An  objective  function  with  the  voxel- 
dependent  prescription  is  constructed  with  incorporation  of  the  nonuniform  dose 
prescription.  The  functional  unit  density  distribution  in  a  sensitive  structure  is  also 
considered  phenomenologically  when  constructing  the  objective  function.  Two  cases 
with  different  hypothetical  biology  distributions  are  used  to  illustrate  the  new  inverse 
planning  formalism.  For  comparison,  treatments  with  a  few  uniform  dose  prescriptions 
and  a  simultaneous  integrated  boost  are  also  planned.  The  biological  indices,  TCP  and 
NTCP,  are  calculated  for  both  types  of  plans  and  the  superiority  of  the  proposed 
technique  over  the  conventional  dose  escalation  scheme  is  demonstrated.  Our 
calculations  revealed  that  it  is  technically  feasible  to  produce  deliberately  nonuniform 
dose  distributions  with  consideration  of  biological  information.  Compared  with  the 
conventional  dose  escalation  schemes,  the  new  technique  is  capable  of  generating 
biologically  conformal  IMRT  plans  that  significantly  improve  the  TCP  while  reducing  or 
keeping  the  NTCPs  at  their  current  levels.  Biologically  conformal  radiation  therapy 
(BCRT)  incorporates  patient  specific  biological  information  and  provides  an  outstanding 
opportunity  for  us  to  truly  individualize  radiation  treatment.  The  proposed  formalism  lays 
a  technical  foundation  for  BCRT  and  allows  us  to  maximally  exploit  the  technical 
capacity  of  IMRT  to  more  intelligently  escalate  the  radiation  dose. 

Key  word:  Inverse  Planning,  Biological  Model,  TCP,  NTCP,  IMRT 
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I.  INTRODUCTION 

Intensity  modulated  radiation  therapy  (IMRT)  has  been  used  clinically  to  provide  a 
highly  conformal  radiation  dose  to  the  target  volume  while  reducing  the  doses  to  the 
surrounding  sensitive  structures!'^.  The  current  IMRT  inverse  planning  is  typically 
aimed  at  producing  a  homogeneous  target  dose  under  the  assumption  of  uniform  biology 
within  the  target  volume.  In  reality,  it  is  known  that  the  spatial  biology  distributions  in 
most  tumors  and  normal  tissues  are  rarely  homogeneous.  To  maximize  the  efficacy  of 
IMRT,  it  is  desirable  to  take  the  inhomogeneous  biological  information  into  account  and 
to  produce  customized  nonuniform  dose  distributions  on  a  patient  specific  basis.  This 
type  of  radiation  treatment  is  referred  to  as  biologically  conformal  radiotherapy 
(BCRT)!4-19  The  simultaneous  integrated  boost  (SIB)  to  elective  volumes  recently 
appeared  in  the  literature!  1,  17,  20  represents  a  simple  example  of  BCRT.  However,  an 
underlying  deficiency  of  the  current  SIB  approach  is  that  the  boost  doses  are  based  on 
previous  experience,  not  patient-specific  biological  information  characterizing  the  spatial 
tumor  burden  distribution. 

To  establish  the  BCRT  treatment  planning  scheme,  three  major  aspects  must  be 
addressed:  (i)  Determination  of  the  distribution  of  biological  properties  of  the  tumor  and 
critical  structures;  (ii)  Prescription  of  the  desired  dose  distribution  for  inverse  planning; 
and  (iii)  Inverse  planning  to  generate  most  faithfully  the  prescribed  nonuniform  dose 
distribution.  Recently  spurred  efforts  in  biological  imaging,  such  as  positron  emission 
tomography  (PET),  single  photon  emission  computed  tomography  (SPECT),  and 
magnetic  resonance  spectroscopy  imaging  (MRSI),  are  aimed  at  providing  solutions  to 
the  first  problem21-31.  To  give  a  few  examples,  the  clonogen  density  in  malignant 
glioma  can  be  obtained  based  on  the  choline/creatine  ratio  through  MRSl29,  30;  tumor 
hypoxia  can  be  quantified  using  PET  imaging  with  fluorinated  misonidazole  (FMISO)27, 
28,  tumor  proliferation  rate  can  be  obtained  based  on  the  voxel  activity  level  in  DNA 
proliferation  imaging  (e.g.,  fluoro-L-thymidine  PET)25,  26,  32?  ancj  iung  functional 

importance  distributions  can  be  obtained  by  perfusion  imaging33.  While  the  development 
of  molecular  imaging  techniques  is  critically  important  in  mapping  out  biology 
distributions,  the  successful  integration  of  this  information  into  IMRT  planning  through 
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steps  (ii)  and  (iii)  is  also  indispensable  to  fully  exploit  the  obtained  biology  information 
to  improve  patient  care.  In  this  study  we  focus  our  efforts  on  the  last  two  problems,  with 
the  optimistic  assumption  that  spatial  biology  distributions  within  a  patient  has  already 
been  determined  from  biological  imaging  or  other  means.  Our  goal  is  to  establish  a 
theoretical  framework  for  quantitatively  incorporating  the  biological  data  into  IMRT 
inverse  optimization  and  to  show  the  advantage  of  the  selective  dose  escalation  scheme  in 
enhancing  tumor  control  probability  (TCP)  and  reducing  the  normal  tissue  complication 
probability  (NTCP).  In  conjunction  with  the  rapid  development  of  molecular  imaging 
techniques,  this  study  lays  a  technical  foundation  for  BCRT  and  provides  a  basis  for 
clinically  realizing  the  new  treatment  strategy  in  the  future. 


II.  METHODS  AND  MATERIALS 

A.  Biological  characterization  and  nonuniform  target  dose  prescription 

We  assume  that  biological  properties  influencing  radiation  treatment  are  characterized 
phenomenologically  by  three  radiobiology  parameters:  clonogen  density  (p), 
radiosensitivity  (a),  and  proliferation  rate  (y).  Generally,  these  parameters  are  voxel 
dependent.  In  this  work  we  concentrate  on  their  spatial  variation  within  tumor,  and  ignore 
the  time  dependence  of  the  last  two  parameters. 

To  accomplish  BCRT,  an  important  step  is  to  derive  the  desired  dose  distribution 
that  maximizes  the  cell  killing  based  on  (p,  a,  y)  metrics.  In  the  case  of  uniform  biology, 
it  is  well  known  that  the  target  dose  should  be  uniformly  distributed.  It  is,  however,  not 
clear  at  all  what  form  of  dose  distribution  should  be  used  to  maximize  the  cell  killing  for 
an  arbitrary  biology  distribution.  We  start  from  a  linear  quadratic  (LQ)  modeP4-36  wjth 
inclusion  of  the  tumor  cell  proliferation.  According  to  this  model,  the  tumor  clonogen 
survival  S,  in  a  voxel  of  volume  F,  after  an  irradiating  dose  D,  is  given  by 
S,  =  p,V,  exp (-«,  £),  +  -AT),  (1) 

where  a ]  =  a,[  1  +  d:  !{ai  / /?,  )] ,  />,  is  the  initial  clonogen  density,  dt  is  the  fractional  dose, 
a,  and  /?,  are  the  linear-quadratic  coefficients  of  the  cell  survival  curve,  yj  =  In  2  IT p  is  the 
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cell  proliferation  rate,  Tp  is  the  potential  cell  doubling  time  and  AT  is  the  overall 
treatment  time.  The  TCP  of  a  voxel  i  can  be  expressed  as 

TCP,  =  exp [-p,V,  exp (-«,£>,  +  y,AT)] .  (2) 

The  TCP  for  the  whole  tumor  is  product  of  the  TCP\  of  all  voxels  within  the  tumor 
volume,  i.e., 

TCP  =  Y\TCp<  •  (3) 

/ 

For  a  given  set  of  {p,  a,  y},  the  task  is  to  find  the  dose  distribution  that  maximizes  the 
TCP.  Because  of  the  limitation  of  normal  tissue  dose  tolerances,  an  arbitrarily  high  dose 
to  the  tumor  cannot  be  achieved  and  certain  constraints  need  to  be  imposed36-41  jn  ijne 
with  previous  researchers3^  37,  40?  we  restrict  the  integral  dose  to  the  tumor  volume  to  a 
constant.  Mathematically,  the  constraint  is  written  as 

5>,A=£,.  (4) 

/ 


where  w,  is  the  mass  of  voxel  i,  and  E,  is  the  integral  target  dose. 

With  the  above  formulation,  the  task  becomes  the  maximization  of  the  TCP  under 
the  constraint  (4).  The  Lagrange  multiplier  method  is  employed  to  solve  the  problem.  In 
this  approach,  a  function 

L(TCPl,...TCPi,...)  =  Y[TCPi+A(Yirr>iDi  -E,)  (5) 


is  introduced,  where  X  is  the  Lagrange  multiplier,  and  the  solution  is  obtained  by  solving 
the  equations 


dL 

8TCPi 


=  0. 


(6) 


When  mass  and  volume  are  equal  for  all  voxels  in  the  target,  using  a  process  similar  to 
Ebert  and  Hoban  40  (see  appendix),  we  obtained  a  general  formula  for  determining  the 
desired  dose,  D‘0  (/') ,  at  the  voxel 


\T  _  <^n’f  _  £) 


D  ( 0  = 


a. 


ref 


1 


a, 


(r, 


ref 


■y^AT  -  A- In 
a 


^  a  ref  P  ref  ^ 


(7) 


v  a'P'  j 

where  Dref  is  the  reference  dose  for  the  voxel  with  reference  radiobiological  parameters 
(pref ,  aref ,  yref  ).  In  general,  Dref  should  be  set  to  a  value  that  yields  a  clinically  sensible 
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TCP  at  the  reference  voxel.  For  a  given  disease  site,  the  radiation  dose  used  in  current 
clinical  practice  with  “intent  to  cure”  can  be  used  as  a  good  starting  point  in  selecting  the 
value  of  Dref.  Using  Eq.  (7),  it  is  straightforward  to  determine  the  desired  target 
prescription  dose  once  the  radiobiological  parameter  ip,  a,  y)  metrics  and  Dref  are  known. 
Note  that  the  desired  dose  distribution  represents  an  ideal  situation  without  considering 
the  specific  dosimetric  tolerances  of  the  sensitive  structures.  In  reality,  this  dose 
distribution  may  or  may  not  be  exactly  realizable.  Nevertheless,  it  sets  a  landmark  and 
serves  as  the  prescription  dose  in  inverse  planning  to  guide  the  dose  optimization  process. 

The  fractional  dose,  d,,  is  required  to  obtain  the  parameter  or,  in  Eq.  (7).  On  the 
other  hand,  dt  is  not  known  until  Z)07  (/)  is  known.  We  use  a  simple  iterative  method  to 
solve  the  dilemma.  First,  the  fractional  dose  is  initially  set  to  dj=Dre/Nf,  Nf  being  the 
fractional  number.  Second,  D0r(z)  is  calculated  using  Eq.  (7)  and  d,=  D^(i)/Nf  is 

updated.  The  new  Dr0  (/)  is  then  obtained  using  the  updated  d,.  We  find  that  D]  (i) 

converges  to  the  solution  in  less  than  5  iterations.  In  this  study  we  set  a/p=10  Gy  for  all 
target  voxels.  The  formalism  proposed  here  is,  however,  general  and  can  be  extended  to 
deal  with  nonuniform  distributions  of  the  a/p  ratio. 

B.  Inverse  planning  with  spatially  nonuniform  dose  prescription 

The  next  logical  step  after  obtaining  the  calculated  prescription  dose  is  to  use  inverse 
planning  to  derive  the  optimal  beam  profiles  that  will  produce  the  prescribed  dose 
distribution.  To  proceed,  we  construct  an  objective  function  to  take  the  known  biological 
information  into  account.  In  addition  to  the  voxel-specific  prescription  as  determined  by 
Eq.  (7),  the  nonlinear  dose  responses  of  tumor  and  normal  structures  are  considered  using 
the  concept  of  equivalent  volume42-48  0f  a  voxel,  which  is  defined  as 
(AFe#),  =^(/) (D(i)/D,)u"  ,  (8) 

where  (A  Veff),  is  the  effective  volume  for  voxel  i  with  volume  V,  and  dose  D(i),  D,  is  the 
desired  dose  for  a  target  voxel  or  the  TD5/5  of  the  corresponding  organ,  and  <j>{i)  is  the 
functional  unit  density.  The  value  of  n  characterizes  the  dose-volume  effect  of  an  organ 
and  reflects  its  architecture  (serial  or  parallel)  of  the  sensitive  structure.  It  is  obtained  by 
fitting  to  clinical  dose-volume  data.  For  a  sensitive  structure,  n  is  a  positive  number  ( n>0 ) 
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while  for  a  target,  n  should  be  assigned  with  a  small  negative  value  (-1  <«<()).  ^(/)  =  1  for 
a  target  voxel. 

A  general  form  of  the  inverse  planning  objective  function  in  the  voxel  domain  is 
written  as 

F = h rr -Trf.  o + iDc  (0  /  Dl  m'"- )  [dc  (o  -  di  o')]2 

r=1  A  r  ,=1 

+  IX  “X  ^  (0[o.  ]’"'•}  A  (02 ,  (9) 

<7=1  I=1 

where  rT  and  ra  are  the  structure  specific  importance  factor  of  target  r  and  sensitive 
structure  a,  respectively,  tT  and  s„  the  number  of  targets  and  sensitive  structures,  Nr  and  Na 
the  total  number  of  voxels  of  target  x  or  sensitive  structure  a,  nT  and  na  the  n  parameter  of 
target  r  and  sensitive  structure  tr,  Dc(i)  the  calculated  dose  in  voxel  i,  Dl  (/)  the 

prescription  dose  in  a  target  voxel  i  given  by  Eq.  (7),  and  TDa^  the  TD5/5  of  sensitive 
structure  a.  The  objective  function  becomes  the  conventional  quadratic  objective  function 
if  the  term  in  the  bracket  inside  each  summation  is  set  to  unity  (this  is  true  when  the  dose- 
volume  effect  is  negligible,  i.e.,  when  nx  -  ns  =+00).  More  detailed  information  about 
the  optimization  algorithm  can  be  found  in  Ref.  49. 

C.  Implementation 

A  software  module  for  optimizing  the  objective  function  (9)  is  implemented  in  the 
platform  of  the  PLUNC  treatment  planning  system  (University  of  North  Carolina,  Chapel 
Hill,  NC).  The  dose  calculation  engine  and  a  variety  of  image/beam/plan  display  and 
evaluation  tools  of  PLUNC  are  used  to  review  and  compare  the  optimization  results.  The 
ray-by-ray  iterative  algorithm  (SIITP)  reported  earlier^*  51  js  employed  to  obtain  the 
optimal  beam  intensity  profiles.  The  dose  volume  histograms  (DVHs)  of  the  involved 
organs  are  displayed  at  the  end  of  each  iterative  step  to  visually  monitor  the  optimization 
process. 

D.  Plan  review  tools 

It  is  desirable  to  extend  the  currently  used  plan  review  tools  to  deal  with  a  biologically 
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heterogeneous  system.  For  a  target,  we  define  the  effective  dose  at  a  voxel  as  the  physical 
dose  normalized  by  the  desired  dose  determined  by  Eq.  (7).  The  effective-dose  volume 
histogram  (EDVH),  which  is  obtained  by  replacing  the  dose  with  the  effective  dose  in 
conventional  DVH,  is  a  useful  tool  for  assessing  BCRT  plans.  For  a  sensitive  structure 
we  replace  the  fractional  volume  by  (Z 5,  V,  to  construct  a  functional  dose  volume  histogram 

(FDVH),  similar  to  that  proposed  by  Lu  et  al^2  and  Marks  et  aP3.  After  including  the 
heterogeneous  biological  information  into  the  EDVH  or  FDVH,  the  wisdom  used  in 
interpreting  a  conventional  DVH  can  be  applied  to  assess  the  BCRT  plans.  In  addition  to 
the  effective  dose  and  the  EDVH  or  FDVH,  a  cluster  of  DVHs,  each  corresponding  to  a 
given  set  of  biological  parameters  {p,  a,  y),  are  also  useful  to  assess  dosimetric  behavior 
of  the  system  as  a  function  of  the  biological  status  of  the  system. 

Besides  the  dosimetric  evaluation  tools,  we  also  used  the  TCP  and  NTCPs  for 
plan  evaluation.  In  calculating  TCP  and  NTCP,  the  heterogeneous  biology  distributions 
need  to  be  taken  into  account.  TCP  is  calculated  using  equations  (2)  and  (3)  and  NTCP  is 
assessed  using  Lyman’s  model.  The  Kutcher-Burman  effective-volume  DVH  reduction 
method44  js  extended  to  include  the  non-uniform  functional  unit  density  distribution 
using  Eq.  (8)  when  transforming  a  nonuniform  dose  distribution  into  a  uniform 

irradiation  of  an  effective  partial  volume.  Model  parameters  from  Burman  et  are 
listed  in  Table  I  for  the  NTCP  calculation. 

E.  Case  studies 

A  prostate  case  with  two  different  hypothetical  distributions  of  radiobiological 
parameters  is  used  to  test  the  proposed  BRCT  inverse  planning  scheme.  In  each  study,  the 
target  consists  of  the  prostate  gland  with  a  few  intra-prostatic  lesions.  The  sensitive 
structures  include  the  rectum,  bladder  and  femoral  heads.  Figures  la  and  3a  show  the 
geometric  shapes  and  locations  of  the  structures  in  the  two  examples. 

In  the  first  example  the  target  includes  four  biologically  different  regions,  and  the 
functional  unit  density  distributions  in  the  sensitive  structures  are  uniform.  Region  1 
represents  the  basis  reference  target  volume  with  typical  parameters^,  55  =  5xl05 

clonogen/cm3,  a,  =0.26  Gy'1,  and  /,  =ln2/40  day’1.  The  radiobiological  parameters  of  the 
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intra-prostatic  lesions  are  listed  in  Table  II.  The  parameters  na  characterizing  the  dose- 
volume  effect  of  the  sensitive  structures  in  the  objective  function  (9)  can  be  found  in 
Table  I.  The  parameter  nr  is  chosen  to  be  -0.2.  For  comparison,  five  IMRT  plans,  indexed 
by  Plan-1,  -2,  -3,  -4,  and  -5,  are  generated.  Plan-1  is  obtained  using  the  BCRT 
optimization  scheme  described  above  with  Dref=  70  Gy.  Plan-2  is  obtained  by  prescribing 
the  whole  target  a  uniform  dose  of  70  Gy.  Plan-3  and  -4  are  similar  to  Plan-2  except  that 

the  dose  is  escalated  to  81  Gy  and  91  Gy  respectively.  Plan-5  is  the  SIB  IMRT 

plan  with  the  same  prescribed  doses  as  that  of  the  BCRT.  In  Plan-1  to  -4,  the  objective 
function  expressed  in  Eq.  (9)  is  used  and  in  Plan-5  the  conventional  dose-based  quadratic 
objective  function  is  adopted.  The  optimization  parameters  (maximum  dose  constraints 
and  importance  factors)  in  the  dose-based  method  were  adjusted  by  trial-and-error  to 
obtain  the  “optimal”  plan.  The  same  beam  configuration  (five  equally  spaced  15MV 
photon  beams  with  gantry  angles  of  0°,  12°,  144°,  216°,  and  288°  in  IEC  convention)  is 
used  in  generating  the  five  plans. 

In  the  second  example  we  hypothetically  introduced  a  higher  functional  unit 
density  region  in  the  rectum  (R_Region  2  as  shown  in  figure  3a)  in  addition  to  three 
biologically  different  target  regions.  The  functional  unit  density  of  the  R  Region  1  is 
assigned  a  value  of  1  and  that  of  the  R_Region  2  is  set  to  be  4.  The  same  set  (p0i ,  a , ,  /,  ) 
as  the  previous  example  and  a  reference  dose  of  70  Gy  are  assigned  to  the  prostate  gland. 
The  parameters  for  other  target  regions  are  listed  in  Table  II.  Once  again,  five  IMRT 
plans  are  generated:  Plan-1  is  obtained  using  the  proposed  selective  dose  escalation 
scheme,  Plan-2,  -3,  and  -4  are  generated  using  different  uniform  prescription  doses  (70, 
81,  and  91  Gy)  and  Plan-5  is  SIB  plan  with  the  same  prescription  as  Plan-1  but  is 
optimized  using  the  conventional  quadratic  objective  function.  In  generating  these  five 
plans,  seven  equally  spaced  15MV  photon  beams  (0°,  51°,  103°,  154°,  206°,  251°,  and 
309°)  are  employed. 


III.  RESULTS 
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A.  Example  1:  Prostate  case  with  four  biologically  different  regions 
In  the  first  example,  based  on  Eq.  (7)  and  the  parameters  listed  in  Table  II,  the 
prescription  doses  to  the  target  region  2,  3  and  4  are  determined  to  be  85  Gy,  1 1 9  Gy  and 
75  Gy,  respectively.  In  order  to  examine  the  capability  of  the  BCRT  inverse  planning 
system  in  producing  an  extremely  nonuniform  dose  within  a  target  volume,  we  have  used 
an  “extreme”  combination  of  {p,  a,  y},  which  leads  to  an  exceedingly  high  prescription 
dose  in  region  3  (119  Gy).  Figures  lb-d  show  the  isodose  distributions  of  Plan-1  in  a 
transverse  slice  and  two  sagittal  slices.  The  EDVH  of  the  target  and  the  DVHs  of  the 
sensitive  structures  are  plotted  in  figure  2  for  Plan-1  in  solid  curves.  For  comparison,  the 
corresponding  EDVHs  and  DVHs  of  Plan-2,  -3,  -4,  and  -5  are  also  shown  in  the  figures 
as  dashed,  dotted,  dash-dotted  and  dash-dot-dotted  curves.  As  seen  from  figure  1,  all 
regions  in  the  prostate  are  well  covered  by  their  prescription  doses  and  the  sensitive 
structures  are  well  spared.  Even  in  this  “extreme”  case,  it  seems  that  the  inverse  planning 
system  can  satisfy  the  biological  requirement.  A  steep  dose  gradient  is  found  at  the 
interface  between  the  target  and  the  rectum.  A  comparison  of  the  target  EDVH  in  figure 
2a  indicates  that  above  98.5%  of  the  target  voxels  achieved  their  desired  doses  in  Plan-1 
and  Plan-5.  However,  for  the  uniform  dose  escalation  scheme,  the  desired  doses  in  some 
regions  (region  2,  3  and  part  of  region  4  in  Plan-2;  region  2  and  3  in  Plan-3;  and  region  3 
in  Plan-4)  are  not  achieved.  We  found  that,  in  Plan-1,  the  doses  to  the  surrounding 
sensitive  structures  are  not  significantly  increased  compared  with  those  of  Plan-2,  despite 
of  the  fact  that  some  voxels  in  region  4  receive  a  dose  as  high  as  1 19  Gy.  In  Plan-1,  the 
rectum,  bladder  and  femoral  heads  are  better  spared  in  comparison  with  Plan-3,  -4. 
However,  by  comparing  the  DVHs  of  Plan-1  and  -5,  it  is  noticed  that,  although  the  target 
coverage  in  Plan-5  is  similar  to  that  in  Plan-1,  the  sensitive  structures  in  Plan-5  receive 
much  higher  doses  than  Plan-1,  indicating  that  the  proposed  approach  can  increase  the 
sensitive  structure  sparing  compared  with  the  conventional  dose-based  quadratic 
objective  function.  In  addition,  as  can  be  expected,  the  target  doses  in  Plan-1  and  -5  are 
less  uniform  in  the  target  volume  in  comparison  with  that  of  Plan-2,  -3,  and  -4.  This  is 
more  pronounced  in  the  target  region  1,  where  about  50%  of  the  volume  receives  a  dose 
larger  than  85  Gy  as  shown  in  figure  2b,  resulting  in  an  effective  doses  above  120%  in 
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-50%  of  the  target  voxels  (see  figure  2a).  However,  the  increase  of  dose  homogeneity  is 
desirable  here  provided  that  the  NTCPs  are  not  compromised. 

Table  III  lists  the  calculated  TCPs  for  the  targets  and  NTCPs  for  the  sensitive 
structures  with  consideration  of  heterogeneous  biology  in  all  plans.  It  is  seen  that  the 
overall  TCPs  for  the  three  plans  with  uniform  target  dose  prescriptions  (Plan-2,  -3,  and  - 
4)  are  all  less  than  that  of  the  BCRT  plan  (Plan-1)  and  SIB  plan  (Plan-5).  This  is 
understandable  because,  in  Plan-2,  -3  and  -4,  some  target  regions  (such  as  target  region 
3)  receive  doses  much  less  than  the  desired  doses.  For  example,  in  Plan-4,  the  TCP  for 
target  region  3  is  only  0.461.  Even  if  the  TCPs  for  region  1, 2,  and  4  are  all  close  to  1.00, 
the  resultant  total  TCP  for  Plan-4  is  0.461.  In  contrast,  the  TCPs  of  Plan-1  and  Plan-5  are 
0.984  and  0.981,  respectively.  Furthermore,  we  found  that  the  NTCPs  of  the  sensitive 
structures  in  Plan-1  are  very  close  to  Plan-2,  significantly  less  than  Plan-3,  -4  and  -5.  For 
example,  the  rectum  NTCPs  are  0.21%  for  Plan-1  and  0.20%  for  Plan-2.  These  are 
increased  to  0.65%,  1.84%  and  0.89%  for  Plan-3,  -4  and  -5,  respectively.  Again,  although 
similar  overall  TCPs  are  achieved  for  the  BCRT  and  dose-based  SIB  IMRT  plans  when 
the  same  dose  prescriptions  are  used,  the  rectum  NTCPs  are  significantly  reduced  using 
the  proposed  formulism.  This  is  consistent  with  our  previous  study  of  the  objective 
function  in  the  context  of  conventional  IMRT  aiming  to  deliver  a  uniform  dose  to  the 
target  volume  (49). 

B.  Example  2:  Prostate  case  with  three  biologically  different  regions  and 
nonuniform  importance  in  rectum 

In  the  second  hypothetical  example,  there  are  three  biologically  different  regions  in  the 
prostate  and  two  functionally  different  regions  in  the  rectum.  The  prescription  doses  for 
the  three  target  regions  are  70  (reference  dose),  99  and  121  Gy,  as  determined  by  Eq.  (7) 
with  the  biological  parameters  listed  in  table  II.  Figures  3b-f  show  the  isodose 
distributions  of  Plan-1  in  three  transverse  slices  and  two  sagittal  slices.  The  EDVHs  and 
DVHs  of  the  target  and  sensitive  structures  for  Plan-1  to  Plan-5  are  plotted  in  figure  4  as 
solid,  dashed,  dotted,  dash-dotted  and  dash-dot-dotted  curves.  Similar  to  the  previous 
example,  in  Plan-1,  all  regions  in  the  prostate  are  well  covered  by  a  dose  comparable  to 
the  prescription  and  the  sensitive  structures  are  well  spared.  The  dose  gradient  at  the 
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interface  between  the  target  and  the  rectum  is  very  sharp  for  all  the  plans.  From  figure  4a 
we  find  that  above  98%  of  the  target  voxels  achieved  their  desired  doses  in  Plan-1 .  As  a 
consequence  of  incorporating  functional  unit  density  information  in  inverse  planning,  the 
rectum  sparing  is  even  better  than  that  of  Plan-2,  much  better  than  that  of  Plan-3,  -4. 
However,  we  notice  that  the  sparing  of  the  femoral  heads  in  Plan-1  is  not  as  good  as  that 
in  Plan-2,  -3  and  -4.  This  is  because  high  intensity  beamlets  that  pass  through  the  femoral 
heads  are  needed  to  adequately  irradiate  the  target  region  3,  as  seen  from  figures  3b  and 
3c.  In  addition,  similar  to  the  first  example,  the  target  coverage  in  Plan-5  is  close  to  that 
in  Plan-1,  but  the  doses  to  the  sensitive  structures  in  Plan-5  are  much  higher  than  that  in 
Plan-1. 

Table  IV  lists  the  calculated  TCPs  and  NTCPs  for  all  plans.  Once  again,  we  found 
that  the  TCP  of  the  target  in  the  proposed  BCRT  technique  is  much  higher  and  the  NTCP 
of  the  rectum  is  lower  compared  with  those  obtained  using  the  conventional  uniform  dose 
escalation  schemes.  Remarkably,  the  overall  TCP  for  the  target  is  increased  from  0.823  to 
0.982  and  the  NTCP  of  the  rectum  is  reduced  from  3.1%  to  0.40%  when  Plan-4  is 
replaced  by  the  selective  dose  escalation  scheme  (Plan-1).  Again,  we  found  that,  for 
similar  overall  TCPs,  the  rectum  NTCPs  of  the  BCRT  plan  are  much  lower  in  comparison 
with  that  obtained  using  dose-based  SIB  scheme. 


IV.  DISCUSSION 


Eq.  (7)  provides  a  general  formula  for  determining  the  desired  target  dose  distribution 
based  on  the  known  biology  information  of  the  system  and  represents  one  of  the  main 
results  of  this  study.  A  few  special  cases  are  worth  discussing  here.  First,  when  the 
biology  distribution  is  uniform  in  the  target,  a  uniform  dose  of  Dref  is  desired.  This  is 
consistent  with  previous  studies^?  and  existing  clinical  knowledge. 

When  the  clonogen  density  p  is  nonuniform  while  the  values  of  a  and  y  are 
constant  across  the  target,  we  have 


Dl{0  =  D«--rln 


(10) 
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which  is  identical  to  the  formula  obtained  by  Webb  and  Nahum36.  Eq.  (10)  indicates  that 
the  desired  dose  depends  on  the  tumor  cell  density  logarithmically  and  is  thus  relatively 
insensitive  to  a  variation  in  the  clonogen  density.  For  <2 't;/.  =0.3 12,  for  example,  even  if  the 

clonogen  density  in  a  tumor  voxel  is  10  times  higher  than  that  of  the  reference  situation, 
the  desired  dose  is  only  about  7  Gy  higher  than  the  reference  value.  A  detailed  discussion 

of  this  special  situation  has  been  presented  by  Webb  and  Nahum36. 

Another  special  case  is  that  the  tumor  clonogen  density  and  the  proliferation  rate 
are  constant  and  the  radiosensitivity  a  is  spatially  nonuniform.  Eq.  (7)  now  becomes 


T  ^ ref  1 

DT  (J)  =  — - 7  In 


The  desired  dose  is  approximately  inversely  proportional  to  the  parameter  a ■  and  is  thus 
sensitively  dependent  on  the  value  of  parameter  a] .  This  is  similar  to  the  conclusions  of 
Ebert  and  Hoban40  and  Levin-Plotnik  and  Hamilton^!.  For  example,  if  a]  is  reduced 
from  0.312  (corresponding  to  a  =0.26,  fractional  dose  d,  =2.0Gy  and  a/p  ratio=10  Gy)  to 
0.18  (corresponding  to  a  =0.15,  <7,  =2.0Gy,  and  a/p  =10  Gy),  the  desired  dose  is  increased 
by  about  70%  (from  70  Gy  to  about  1 1 8  Gy). 

If  we  keep  the  tumor  clonogen  density  and  radiosensitivety  a’  constant  and  only 
allow  the  proliferation  rate  y  to  vary  spatially,  then 

flr(0  =  A*+Ao',->w)A7\  (12) 

Thus  the  desired  dose  increases  linearly  with  the  proliferation  rate.  In  this  work  the 
potential  cell-doubling  times,  Tp,  used  by  King  et  afi 4  are  adopted.  Since  Tp  for  a  prostate 
tumor  is  relatively  longer,  its  influence  on  the  desired  dose  is  not  very  significant. 
However,  for  other  more  rapidly  proliferating  tumors,  the  proliferation  rate  may  play  an 
important  role.  In  such  situations,  reducing  the  overall  treatment  time  AT  (e.g.,  using  an 
accelerated  scheme)  is  helpful  to  minimize  the  influence  of  the  proliferation  rate. 

We  emphasize  that  the  quadratic  term  in  the  linear-quadratic  model  plays  an 
important  role  in  accounting  for  the  fractionation  effect.  If  only  the  linear  term  is  kept,  the 
total  dose  Z>7  (/)  in  Eq.  (7)  is  no  longer  entangled  with  the  fractional  dose  dt.  When  the 
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quadratic  term  is  “switched  on”,  the  value  of  D 7  (/)  depends  not  only  the  total  reference 

dose  but  also  the  fractional  dose.  For  a  large  fractional  dose,  the  total  dose  will  be  less, 
and  vice  versa.  In  other  words,  the  total  dose  received  by  a  voxel  is  determined  by  two 
contributing  factors,  one  being  the  local  biological  parameters  {p,a,y},  and  the  other  being 
the  coupling  between  the  fractional  dose  and  the  total  dose.  The  latter  is  responsible  for 
the  phenomenon  that  the  total  dose  needs  to  be  decreased  when  the  number  of  fractions  is 
reduced.  If  the  quadratic  term  were  ignored,  according  to  Eq.  (7),  the  dose  required  at  a 
voxel  would  be  much  higher.  For  example,  the  desired  doses  for  target  region  3  in 
Example  1  are  determined  to  be  119  Gy  and  135  Gy  with  and  without  inclusion  of  the 
quadratic  term,  respectively. 

We  also  would  like  to  emphasize  that  in  this  study,  the  radiosensitivity  a'  and 
proliferation  rate  y  are  assumed  to  be  constants  during  the  whole  treatment  course.  In 
reality,  both  a'  and  y  may  change  with  time  due  to  such  biological  processes  as  tumor  cell 
redistribution5^  and  reoxygenation5?.  The  time  dependence  of  these  factors  may  result  in 
a  reduction  of  the  desired  prescription  dose,  and  this  effect  should  be  investigated  in  the 
future. 

Comparing  with  the  uniform  dose  escalation  scheme,  our  study  clearly  suggests  that 
deliberately  incorporating  an  inhomogeneous  dose  distribution  significantly  enhances  the  TCP 
and  reduces  the  NTCP.  Physically,  we  believe  that  the  significant  improvement  arises  from 
the  more  effective  use  of  radiation  in  the  newly  proposed  treatment  scheme.  A  great  deal  of 
dose  is  “wasted”  in  the  conventional  uniform  dose  escalation  scheme.  For  example,  in  the  first 
example  the  increased  doses  in  the  target  region  1  and  4  have  almost  no  contributions  to  the 
enhancement  of  the  TCP  when  Plan-2  (70  Gy  uniform  dose  to  the  prostate  gland)  is  replaced 
by  Plan-3  (81  Gy)  or  Plan-4  (91  Gy).  Even  though  part  of  the  prostate  receives  high  doses  in 
the  selective  dose  escalation  scheme  (for  example,  119  Gy  in  target  region  3  of  the  first 
example),  the  total  deposited  energy  in  the  targets  is  still  less  than  that  of  Plan-3  or  -4.  It  is 
thus  not  difficult  to  understand  why  deliberately  nonuniform  dose  distributions  can,  in  general, 
reduce  the  radiation  side-effects  and  represent  a  more  intelligent  way  to  irradiate  the  tumor 
target. 

A  similar  deficiency  also  exists  in  the  current  SIB  approach.  Although  it  is  clear  that 
the  regions  with  different  tumor  burdens  should  be  given  different  doses,  the  specific 
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values  for  the  regions  are  determined  in  an  ad  hoc  manner.  The  empirical  boost  dose  could 
be  too  low,  in  which  case  the  tumor  control  is  sacrificed,  or  too  high,  in  which  case  other 
parts  of  the  system  are  compromised.  The  problem  is  aggravated  when  the  tumor  burden 
varies  continuously  from  point  to  point.  In  the  proposed  BCRT  approach,  the  prescribed 
dose  is  voxel-dependent  and  is  determined  based  on  the  tumor  biology  distribution.  In 
addition,  a  more  sophisticated  objective  function  is  developed  to  take  the  dose  volume 
effect  and  functional  density  information  of  the  sensitive  structures  into  account,  resulting 
in  better  sparing  of  the  sensitive  structures. 

Finally,  it  should  be  recognized  that  our  knowledge  of  radiobiological  parameters 
for  tumors  or  normal  tissues  is  still  very  crude  and  the  validity  of  the  model  is  still  under 
establishment.  Therefore,  the  LQ  model  and  the  parameters  adopted  in  the  paper  are  fine 
for  a  proof  of  principle  but  they  should  not  be  taken  as  more  than  that. 


V.  CONCLUSION 

In  the  presence  of  nonuniform  biology  distributions,  IMRT  inverse  planning  is  complicated  by 
the  fact  that  it  is  not  clear  what  represents  the  appropriate  spatial  dose  prescription,  which  is 
generally  used  as  a  landmark  to  guide  the  dose  optimization  process.  In  this  work,  we  have 
described  a  technique  for  deriving  the  prescription  dose  based  on  a  LQ  model  with 
consideration  of  the  cell  proliferation.  The  relation  is  quite  general  and  can  be  used  as 
prescription  dose  to  guide  an  arbitrary  inverse  planning  objective  function  aimed  to  produce 
customized  dose  distribution  in  accordance  with  the  spatial  biology  information.  For  a  given 
patient,  IMRT  inverse  planning  now  consists  of  two  steps:  Derivation  of  the  prescription  dose, 
and  beam  profile  optimization  that  produces  as  closely  as  possible  this  prescription  dose.  The 
formalism  proposed  here  lays  a  technical  foundation  for  future  BCRT  development,  allowing 
us  to  escalate  tumor  dose  more  intelligently  and  effectively.  When  combined  with  state-of-the- 
art  biological  imaging  techniques,  which  promise  to  reveal  detailed  patient-specific  biology 
distribution  information,  this  study  may  have  significant  implication  for  the  management  of 
cancer  in  the  future. 
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Appendix 

We  present  the  detailed  derivation  process  for  Eq.  (7)  under  the  condition  of  equal  mass 
and  volume  for  all  target  voxels. 

Substituting  Eq.  (5)  into  the  Eq.  (6),  we  obtain, 


XTCP.  =  ATCP.„,  =  _TCp  _ 


'  8TCP 


ref 


8TCP. 


(Al) 


Since  X  ^  0,  otherwise,  TCP  becomes  zero  according  to  the  requirement  of  Eq.  (6), 
which  corresponds  to  a  minimum.  If  we  assumed  that  mass  for  all  target  voxels  is  equal, 
then  Eq.  (Al)  becomes 


TCP.  -  dD'  -  =  TCP. 


8D 


ref 


ref 


8TCPj 
From  Eq.  (2)  we  have 

Dt  =  — InJ  —  /, AT  +  In 

a.  | 


8TCP. 


(A2) 


-1 

pXi 


ref 


In  TCP 


(A3) 


Substituting  the  expressions  from  Eq.  (A3)  for  both  D,  and  A*/ into  Eq.  (A2),  we  have 
a,\n(TCp)  =  aref\n(TCPref).  (A4) 

The  desired  doses,  D ]  (/) ,  producing  maximum  TCP  with  the  constraint  of  constant 

integral  dose,  can  be  obtained  by  substituting  TCP,  and  TCPref  expressed  in  Eq.  (2)  into 
Eq.  (A4) 


Z>; (i)  =~~~D  f  -r,)A7- In 

a. 


a, 


a. 


aref  Pref^ref 

aiPlVi 


(A5) 


When  volume  for  all  target  voxels  is  equal,  Eq.  (A5)  becomes  Eq.  (7). 
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Table  I  Dose-volume  parameters  of  various  sensitive  structures  used  for  calculating 
NTCP  in  this  study. 


Sensitive  structures 

n 

m 

Dsns  (Gy) 

Bladder 

0.50 

0.11 

80 

Rectum 

0.12 

0.15 

80 

Femoral  head 

0.25 

0.12 

65 

Table  II  Radiological  parameters  for  the  target  regions  in  the  two  examples. 


Targets 

p0j  (clonogen/cm3) 

ai  (Gy ') 

Yi(  day1) 

Region  2 

~  5x10s 

0.26 

ln2/40 

Example  1 

Region  3 

5xl05 

0.13 

ln2/40 

Region  4 

5x10s 

0.26 

ln2/10 

Example  2 

Region  2 
Region  3 

5xl06 

5xl03 

0.20 

0.10 

■nv 
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Table  III  Comparison  of  TCP  and  NTCP  for  the  four  IMRT  plans  for  Example  1 . 


Plan-1 

Plan-2 

Plan-3 

Plan-4 

Plan-5 

(BCRT  plan)  (70  Gy  Uniform) 

(81  Gy  Uniform) 

(91  Gy  Uniform) 

(SIB  plan) 

Region  1 

0.997 

0.995 

1.000 

1.000 

0.994 

TCP 

Region  2 

0.998 

0.642 

0.995 

1.000 

0.999 

Region  3 

0.989 

0.000 

0.002 

0.461 

0.989 

Region  4 

1.000 

0.997 

1.000 

1.000 

0.998 

Overall 

0.984 

0.000 

0.002 

0.461 

0.981 

Rectum 

0.212 

0.196 

0.652 

1.84 

0.885 

NTCP 

Bladder 

1.6x1 0'5 

1.4x1  O'5 

2.3x1  O'5 

4.2x1 0’5 

3.6x1 0‘5 

(%) 

Femoral  head  (R) 

2.0x1 0‘5 

2.1  xlO'6 

2.6x1 O'4 

1.75xl0"4 

3.9x1 0‘5 

Femoral  head  (L) 

1.2x1 0'5 

2.0x1 0‘6 

7.0x1 0‘4 

5.26x1 0‘4 

6.9x1 0'5 

Table  IV  Comparison  of  TCP  and  NTCP  for  the  four  IMRT  plans  for  Example  2. 


Plan-1 

(BCRT  plan) 

Plan-2 

(70  Gy  Uniform) 

Plan-3 

(81  Gy  Uniform) 

Plan-4 

(91  Gy  Uniform) 

Plan-5 

(SIB  Plan) 

Region  1 

0.997 

0.995 

1.000 

1.000 

0.968 

TCP 

Region  2 

0.989 

0.000 

0.587 

0.981 

0.987 

Region  3 

0.996 

0.006 

0.408 

0.839 

0.990 

Overall 

0.981 

0.000 

0.239 

0.823 

0.946 

Rectum 

0.397 

0.414 

1.46 

3.12 

1.25 

NTCP 

Bladder 

1.5xl0'5 

1.2x1 0'5 

1.8x1  O'5 

4.3x1  O'5 

3.9x1  O'5 

(%) 

Femoral  head  (R) 

3.7X  1 0"5 

1.5x1  O'5 

1.8x1  O'5 

5.3x1 0'5 

2.3x1  O'5 

Femoral  head  (L) 

4.9x1 0'5 

1 . 1  x  1  O’5 

3.0x1 0'5 

4.5x1  O’5 

3.6x1  O'5 
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Legends 


Figure  1.  A  hypothetical  prostate  case  with  four  biologically  different  regions  (Example  1).  (a) 
Geometric  shapes  and  locations  of  the  targets  and  sensitive  structures;  (b)-(d):  Isodose 
distributions  in  an  axial  slice  and  two  sagittal  slices  for  Plan-1,  generated  by  optimizing  the 
objective  function  with  a  nonuniform  dose  prescription  derived  from  Eq.  (7). 

Figure  2.  Comparison  of  EDVHs  and  DVHs  of  the  BCRT  plan  (Plan-1),  three  uniform  1MRT 
plans  (Plan-2:  70Gy,  Plan-3:  81  Gy,  and  Plan-4:  91  Gy)  and  the  SIB  plan  (Plan-5)  in  example  1. 

(a) :  Target  EDVHs  for  the  five  plans  (Insert  is  the  regular  DVHs  of  the  prostate  target).  The 
normalized  doses  to  the  target  region  1,  2,  3  and  4  are  70,  85,  119  and  75  Gy,  respectively;  (b)- 
(e):  DVHs  of  different  target  regions  and  sensitive  structures  for  the  five  plans.  The  solid,  dashed, 
dotted,  dash-dotted  and  dash-dot-dotted  curves  represent  the  results  of  Plan- 1,  -2,  -3,  -4  and  -5, 
respectively.  The  effective  dose  is  defined  as  the  physical  dose  at  a  voxel  normalized  by  its 
desired  dose  determined  by  Eq.  (7). 

Figure  3.  A  hypothetical  prostate  case  with  three  biologically  different  regions  and  nonuniform 
importance  in  the  rectum  (Example  2).  (a)  Geometric  shapes  and  locations  of  the  targets  and 
sensitive  structures;  (b)-(d):  Isodose  distributions  in  three  axial  slices  and  two  sagittal  slices  for 
Plan-1,  generated  by  optimizing  the  objective  function  with  nonuniform  dose  prescription  derived 
from  Eq.  (7). 

Figure  4.  Comparison  of  EDVHs,  FDVHs  and  DVHs  of  the  BCRT  plan  (Plan-1),  three  uniform 
IMRT  plans  (Plan-2:  70Gy,  Plan-3:  81Gy,  and  Plan-4:  91Gy)  and  the  SIB  plan  (Plan-5)  in 
example  2.  (a):The  target  EDVHs  for  the  five  plans  (Insert  is  the  regular  DVHs  of  the  prostate 
target).  The  normalized  doses  to  the  target  region  1,  2,  and  3  are  70,  99  and  121  Gy,  respectively; 

(b) :  The  rectum  FDVHs  for  the  five  plans  (Insert  is  the  regular  DVHs  of  the  rectum);  (c)  -(e): 
DVHs  of  different  target  regions  and  sensitive  structures  for  the  five  plans.  The  solid,  dashed, 
dotted,  dash-dotted  and  dash-dot-dotted  curves  represent  the  results  of  Plan-1,  -2,  -3,  -4  and  -5, 
respectively. 
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Abstract 

Clinical  IMRT  treatment  plans  are  currently  made  using  dose-based 
optimization  algorithms,  which  do  not  consider  the  nonlinear  dose-volume 
effects  for  tumours  and  normal  structures.  The  choice  of  structure  specific 
importance  factors  represents  an  additional  degree  of  freedom  of  the  system 
and  makes  rigorous  optimization  intractable.  The  purpose  of  this  work  is  to 
circumvent  the  two  problems  by  developing  a  biologically  more  sensible  yet 
clinically  practical  inverse  planning  framework.  To  implement  this,  the  dose- 
volume  status  of  a  structure  was  characterized  by  using  the  effective  volume 
in  the  voxel  domain.  A  new  objective  function  was  constructed  with  the 
incorporation  of  the  volumetric  information  of  the  system  so  that  the  figure  of 
merit  of  a  given  IMRT  plan  depends  not  only  on  the  dose  deviation  from  the 
desired  distribution  but  also  the  dose-volume  status  of  the  involved  organs.  The 
conventional  importance  factor  of  an  organ  was  written  into  a  product  of  two 
components:  (i)  a  generic  importance  that  parametrizes  the  relative  importance 
of  the  organs  in  the  ideal  situation  when  the  goals  for  all  the  organs  are  met; 
(ii)  a  dose-dependent  factor  that  quantifies  our  level  of  clinical/dosimetric 
satisfaction  for  a  given  plan.  The  generic  importance  can  be  determined 
a  priori,  and  in  most  circumstances,  does  not  need  adjustment,  whereas  the 
second  one,  which  is  responsible  for  the  intractable  behaviour  of  the  trade¬ 
off  seen  in  conventional  inverse  planning,  was  determined  automatically.  An 
inverse  planning  module  based  on  the  proposed  formalism  was  implemented 
and  applied  to  a  prostate  case  and  a  head-neck  case.  A  comparison  with  the 
conventional  inverse  planning  technique  indicated  that,  for  the  same  target 
dose  coverage,  the  critical  structure  sparing  was  substantially  improved  for 
both  cases.  The  incorporation  of  clinical  knowledge  allows  us  to  obtain 
better  IMRT  plans  and  makes  it  possible  to  auto-select  the  importance  factors, 
greatly  facilitating  the  inverse  planning  process.  The  new  formalism  proposed 
also  reveals  the  relationship  between  different  inverse  planning  schemes  and 
gives  important  insight  into  the  problem  of  therapeutic  plan  optimization. 


003 1-9155/04/225 101+17S30. 00  ©  2004  IOP  Publishing  Ltd  Printed  in  the  UK 


5101 


5102 


Y  Yang  and  L  Xing 


In  particular,  we  show  that  the  EUD-based  optimization  is  a  special  case  of  the 
general  inverse  planning  formalism  described  in  this  paper. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version) 


1.  Introduction 

An  important  issue  in  inverse  planning  is  how  to  formalize  the  clinical  goals  to  objectively 
evaluate  the  figures  of  merit  of  different  IMRT  plans.  Despite  intense  research  effort  in 
modelling  the  clinical  decision-making  strategies  (Amols  and  Ling  2002,  Deasy  et  al  2002, 
Earl  et  al  2003,  Hou  et  al  2003,  Lahanas  et  al  2003,  Langer  et  al  1993,  1998,  Lee  et  al 
2000,  Llacer  et  al  2001,  Mohan  et  al  1994,  Webb  2004,  Xing  et  al  1999,  Yan  et  al  2003), 
the  appropriate  form  of  the  objective  function  remains  illusive.  Presently,  two  types  of 
objective  functions  are  widely  used:  dose  or  dose-volume  histogram  (DVH)-based  (physical 
objective  functions)  (Chen  et  al  2002,  Cho  et  al  1998,  Holmes  et  al  1995,  Hristov  et  al 
2002,  Michalski  et  al  2004,  Starkschall  et  a!  2001,  Shepard  et  al  2002,  Xing  et  al  1998) 
and  dose-response-based  objective  functions  (biological  objective  functions)  (Brahme  2001, 
Kallman  et  al  1992,  Miften  et  al  2004,  Mohan  et  al  1 992,  Wang  et  al  1 995,  Webb  and  Nahum 
1993).  The  underlying  difference  between  these  models  lies  in  what  endpoints  are  used  to 
evaluate  the  treatment  or  which  fundamental  quantities  are  used  to  define  the  optimality.  The 
physical  approach  emphasizes  the  difference  between  the  calculated  and  prescribed  doses  and 
does  not  consider  the  nonlinear  effects  for  tumours  and  normal  structures.  Dose-volume 
constraints  are  often  introduced  to  select  a  solution  with  certain  shapes  of  the  DVHs  for  the 
target  and  sensitive  structures.  However,  it  is  important  to  note  that  the  construction  of  the 
DVH  constraints  is  a  priori  in  nature.  The  use  of  constraints  can  only  passively  restrict 
the  accessible  DVHs  by  narrowing  the  solution  space,  and  the  figures  of  merits  of  the  physically 
realizable  plans  are  not  changed  as  long  as  they  satisfy  the  constraints.  To  reflect  our  preference 
over  certain  DVHs  for  a  structure,  it  is  necessary  to  express  the  objective  of  the  structure  as  a 
function  of  the  volumetric  status,  which  has  not  been  achieved  up  to  this  point.  On  the  other 
hand,  biological  model-based  optimization  proponents  argue  that  plan  optimization  should 
be  guided  by  estimates  of  biological  effects,  which  depend  on  the  radiation  dose  through  the 
dose-response  function.  The  treatment  objective  in  biological  model-based  inverse  planning  is 
usually  stated  as  the  maximization  of  the  tumour  control  probability  (TCP)  while  maintaining 
the  normal  tissue  complication  probability  (NTCP)  to  within  acceptable  levels  (Brahme  2001, 
Kallman  et  al  1992,  Langer  et  al  1998,  Mohan  et  al  1992,  Wang  et  al  1995).  In  principle, 
the  biologically  based  models  are  most  relevant  for  radiotherapy  plan  ranking.  However,  the 
dose-response  function  of  various  structures  is  not  sufficiently  understood  and,  at  this  point, 
there  is  considerable  controversy  about  the  models  for  computing  dose-response  indices 
and  their  use  in  optimization.  In  reality,  the  use  of  dose-response  indices  for  optimization 
may  lead  to  very  inhomogeneous  target  dose  distributions.  Furthermore,  it  is  difficult  for 
clinicians  to  specify  the  optimization  criteria  in  terms  of  dose-response  indices.  This  becomes 
compounded  when  two  or  more  independently  optimized  plans  are  to  be  combined.  Because 
of  the  problems,  the  use  of  biological  model-based  optimization  has  mainly  been  limited  in 
the  research  community  and  little  effort  has  been  made  to  implement  the  model  in  commercial 
planning  systems.  Given  the  fact  that  biological  outcome  is  the  ultimate  endpoint  of  radiation 
therapy,  the  importance  of  the  biological  modelling  and  the  biological  model-based  inverse 
planning  can  never  be  overstated. 
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To  pin  down  the  underlying  problem  of  the  current  inverse  planning  formalism  and 
illustrate  the  need  for  a  clinically  more  relevant  approach,  let  us  take  parotid  glands  as  an 
example.  It  is  well  known  that  the  clinical  endpoint  is  the  same  if  the  glands  are  irradiated 
with  15  Gy  to  67%,  or  30  Gy  to  45%,  or  45  Gy  to  24%  of  the  total  volume  (Eisbruch  etal 
1999).  If  a  dose-based  metric,  such  as  the  commonly  used  quadratic  objective  function,  is 
used,  the  rankings  for  the  three  different  scenarios  would  be  completely  different.  Even  with 
the  use  of  dose-volume  constraints,  it  is  difficult,  if  not  impossible,  to  incorporate  this  type 
of  knowledge  to  correctly  model  the  behaviour  of  the  organ  in  response  to  radiation.  Indeed, 
a  constraint  in  optimization  acts  as  a  ‘boundary  condition’  during  the  optimization  and  does 
not  change  the  rankings  of  dosimetrically  different  plans.  This  example  clearly  reveals  the 
inadequacy  of  the  conventional  dose-based  objective  function  and  suggests  the  urgent  need 
for  a  clinically  more  sensible  model.  Obviously,  a  minimum  requirement  for  the  model  is 
that  it  should  be  consistent  with  the  existing  clinical  outcome  data.  For  parotid  glands,  for 
instance,  the  three  different  DVHs  mentioned  above  should  be  scored  equally.  This  type  of 
‘degeneracy’  can  be  achieved  by  effectively  integrating  clinical  endpoint  data  into  the  inverse 
planning  formulation.  For  a  given  patient,  the  specific  DVH  selection  will  be  determined  by 
the  optimization  algorithm  with  the  consideration  of  the  dosimetric/clinical  requirements  of 
other  structures  and  the  results  will,  of  course,  depend  on  the  geometric  and  dosimetric  details 
of  the  given  patient.  The  example,  however,  underscores  the  important  role  of  the  existing 
clinical  data  in  inverse  planning  and  emphasizes  the  essential  ingredients  for  a  clinically 
realistic  objective  function.  Towards  developing  a  biologically  more  sensible  yet  clinically 
practical  inverse  planning  formalism,  in  the  following  we  propose  a  method  to  incorporate 
clinical  endpoint  data  into  the  construction  of  the  objective  function  and  attempt  to  bridge 
the  gap  between  the  clinical  decision-making  process  and  the  computational  modelling.  Our 
study  indicates  that  the  clinical  knowledge-based  modelling  allows  us  to  objectively  rank 
IMRT  plans  according  to  their  clinical  merits  and  makes  it  possible  to  obtain  truly  optimal 
IMRT  plans  with  much  reduced  efforts. 


2.  Methods  and  materials 

2.1.  Dose-volume  effect 

Generally,  the  dose  response  of  a  structure  with  respect  to  the  irradiated  dose  and  volume  is 
complicated.  The  fact  that  the  dose  distribution  in  tumour  or  a  sensitive  structure  is  generally 
inhomogeneous  makes  the  establishment  of  such  a  relationship  even  more  intractable.  Over  the 
last  two  decades,  attempts  have  been  made  by  many  researchers  to  capture  the  main  feature(s) 
of  the  dose-volume  effects.  A  power-law  model  represents  one  of  the  successful  techniques 
in  dealing  with  the  dose-volume  effects  of  sensitive  structures  (Lyman  and  Wolbarst  1987). 
In  this  model  an  equivalent  dose  uniformly  irradiating  the  whole  organ,  Dcq,  can  be  used 
to  represent  the  situation  in  which  a  fractional  partial  volume,  v,  is  irradiated  to  a  dose,  D, 
by  a  simple  power-law  model:  Dcq  =  vll"D.  A  remarkable  characteristic  of  this  model  is 
that,  although  only  a  single  organ-specific  parameter,  n,  is  used,  clinical  and  biological  data 
have  shown  that  this  power  law  holds  well  at  low  complication  levels  (Lyman  and  Wolbarst 
1987,  Schultheiss  etal  1983).  On  the  basis  of  this  relation,  Mohan  etal  (1992)  introduced  the 
concept  of  effective  dose  to  represent  a  non-uniform  dose  distribution  in  a  sensitive  structure. 
Kutcher  and  Burman  (1989)  applied  the  same  power  model  independently  to  each  volume 
element  of  the  histogram  and  introduced  the  concept  of  effective  volume  to  reduce  the  DVH 
of  an  inhomogeneous  dose  distribution  in  a  sensitive  structure  to  a  uniform  dose  distribution. 
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Following  their  study,  in  this  work  we  define  the  effective  volume  (A  VCff);  for  a  voxel  i  with 
volume  A  V  and  dose  Dc(i)  as  follows: 

(AVcff);  =  AV  (Dc{i)/  Acf)'7”  (1) 

and  extend  this  concept  to  handle  the  voxels  in  the  tumour  target,  where  n  is  an  organ-dependent 
parameter,  and  Drcf  is  the  reference  dose.  For  a  sensitive  structure,  n  is  a  small  positive  number 
(0  <  n  <  1)  and  the  value  of  parameter  n  reflects  the  architecture  (serial  or  parallel)  of  the 
sensitive  structure.  For  a  target,  n  should  be  assigned  with  a  small  negative  value  (—  1  < 
n  <  0).  The  biological  meaning  of  equation  (1)  is  that  for  a  sensitive  structure  a  small  volume 
receiving  a  higher  dose  than  a  reference  dose  would  be  equivalent  to  a  larger  volume  receiving 
the  reference  dose;  for  a  target,  a  volume  with  a  lower  dose  would  have  a  larger  effective 
volume.  The  effective  volumes  of  all  voxels  reflect  the  DVH  status  of  the  given  organ,  and  for 
inverse  planning,  this  permits  us  to  deal  with  the  complicated  dose-volume  effect  in  the  voxel 
domain. 


2.2.  Dose-volume-based  objective  function 


The  objective  function,  /,  expressed  as  a  function  of  the  effective  volume  in  the  voxel  domain 
for  an  organ  should  take  the  form  of 

/  =  /({(AVcff),}).  (2) 


Generally,  the  dose-volume  effect  suggests  that  the  voxels  receiving  different  doses  are 
inequivalent:  the  one  with  a  larger  effective  volume  (higher  dose  for  a  critical  organ)  should 
be  penalized  more  when  compared  to  a  voxel  with  a  smaller  effective  volume  (lower  dose). 
Thus  we  heuristically  write  the/  in  the  following  form: 


/  =  1  +  mYlr‘[Dc(i)/D'c(]'/n  +  m  r>>[D[.(/)/Drcf]1/n 

i  l  i 


(3) 


where  equation  (1)  has  been  incorporated,  r,  is  the  importance  factor  of  the  ;'th  voxel, 
representing  the  intrastructural  trade-off  due  to  physical/clinical  requirements  other  than 
the  dose-volume-based  penalty,  /ji  and  are  phenomenological  parameters  of  the  model.  In 
equation  (3),  the  third  (and  higher  order)  term  emphasizes  more  the  voxels  with  high  effective 
volumes,  whereas  the  first  and  second  terms  ensure  that  the  voxels  with  low  effective  volumes 
receive  an  adequate  penalty.  We  typically  set  r,-  s  1,  unless  there  are  other  physical /clinical 
considerations  (e.g.,  when  the  density  of  clonogens  varies  spatially  (Xing  et  al  2002)).  In  this 
work,  we  set  ?)i  =  1  and  =  173  =  •  ■  •  =  0. 


2.3.  Hybrid  of  dose-based  and  dose-volume-based  objective  functions 

Equation  (3)  provides  a  good  description  of  the  dose-volume  effect.  With  proper  choice  of  the 
parameter,  n,  the  clinically  observed  dose-volume  effect  can  be  reproduced  by  the  objective 
function.  In  reality,  other  requirements,  such  as  the  target  dose  uniformity,  should  also  be 
considered.  A  more  general  form  of  inverse  planning  objective  function  can  be  written  as  a 
hybrid  of  the  dose-volume-based  and  the  dose-based  functions.  In  this  situation,  the  overall 
objective  function  of  the  system  takes  the  form  of 
r,  ,  n, 

F  =  {1  +i[Dc«)/DT<Tcf]l'n'}\Dc(i)  -  D0r(i) |*' 

T— 1  T  i=  I 

■fa  1  Mr 

+  E{1-<[Of(i)/Dff,rcf]l/-}Dc(/A, 

<7  =  1  °  /=! 


(4) 
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where  tr  and  s„  are  the  numbers  of  targets  and  sensitive  structures,  Dj  (i)  is  the  prescription 
dose  in  target  voxel  subscripts  r  and  a  represent  target  r  and  sensitive  structure  a,  Nt, 
N„,  rx,  ra<  nx,  n„,  Dr  rcf,  D„Jcf,  kx  and  k„  represent  the  total  numbers  of  voxels,  structure 
specific  importance  factors,  n  parameters,  reference  doses,  power  of  dosimetric  deviation 
from  the  specified  criteria  for  target  z  and  sensitive  structure  a,  respectively.  The  factor 
\Dc(i)  —  Dq  (i')f'  for  target  or  Dc(i)k°  for  a  sensitive  structure  represents  the  contribution 
from  dosimetric  deviation  from  the  ideal  situation.  If  the  kx  and  ka  are  set  to  zero,  the  objective 
function  becomes  purely  dose-volume  driven.  In  particular,  if  we  set  ka  to  zero  and  kx  to  a 
nonzero  value,  the  objective  function  for  a  target  becomes  a  hybrid  of  dose-volume  and  dose, 
whereas  the  objective  functions  for  critical  structures  remain  purely  dose-volume  based.  On 
the  other  hand,  when  all  the  n  parameters  in  equation  (4)  are  set  to  be  +oo,  no  dose-volume 
effects  are  considered  and  equation  (4)  is  reduced  to  the  conventional  dose-based  objective 
function. 


2.4.  Automatic  determination  of  structure  specific  importance  factors 

The  selection  of  structure  specific  importance  factors,  rr  or  r a  in  equation  (4),  is  generally 
done  empirically  by  trial  and  error.  Here  we  describe  an  automated  approach  for  solving 
the  problem.  The  key  to  success  is  to  establish  an  effective  method  to  express  the  structural 
importance  factor  in  terms  of  physically  or  clinically  more  meaningful  quantities.  For  this 
purpose,  we  write  the  importance  of  a  structure  into  a  product  of  two  components:  (i)  a 
generic  factor  parametrizing  the  relative  importance  of  the  organs  in  an  ideal  situation  when 
the  goals  for  the  organs  are  met;  and  (ii)  a  dose-dependent  factor  quantifying  our  level 
of  clinical/dosimetric  satisfaction  for  a  given  plan.  The  first  factor  can  be  determined 
a  priori,  and  in  most  circumstances,  does  not  need  adjustment  (generally  speaking,  the 
value  of  r*  is  determined  on  the  basis  of  the  treatment  modality  and  the  patient’s  overall 
condition),  whereas  the  second  one  is  responsible  for  the  intractable  behaviour  of  the  trade¬ 
off  in  conventional  planning  and  can  be  automatically  determined.  This'  decomposition  is 
essentially  to  normalize  the  conventional  importance  factor  in  terms  of  our  clinical  goals  for 
the  structures  under  discussion.  Because  of  this  decomposition,  the  meaning  of  the  importance 
factor  becomes  more  transparent  and  the  determination  of  the  factors  becomes  straightforward. 
Mathematically,  we  write  ra  =  where  represents  the  first  contribution  described  above 
(the  desired  weighting  among  different  structures  in  an  ideal  situation),  and  r *  is  the  second 
component  and  is  defined  as  a  function  of  NTCP  for  a  sensitive  structure.  In  this  study  r* 
was  set  empirically  (see  tables  2  and  3  for  examples).  r°  is  updated  according  to  the  DVH 
or  the  dose  distribution  during  the  optimization  process  and  reflects  the  most  current  status  of 
trade-off  in  the  system.  Generally,  the  importance  of  a  sensitive  structure  should  be  increased 
in  the  next  iteration  if  NTCP  is  high,  and  vice  versa.  We  found  that  a  simple  linear  relation 
between  and  NTCP, 

ra  =  NTCPff  +  8,  (5) 

describes  the  trade-off  behaviour  of  the  system  well,  where  the  value  of  NTCPff  depends  on 
the  dose  distribution  at  the  current  iteration  for  structure  a ,  and  8  is  a  cut-off  factor  for  NTCP, 
which  is  introduced  to  ensure  the  sensitive  structure  receives  a  minimum  penalty  even  if  its 
NTCP  is  close  to  zero.  We  set  8  as  an  organ-independent  constant  of  0.01%. 

The  NTCP  was  assessed  using  Lyman’s  model  in  this  study.  For  non-uniform  irradiation, 
the  Kutcher-Burman  (1989)  effective-volume  DVH  reduction  method  is  used  to  transform 
a  DVH  into  a  uniform  irradiation  on  an  effective  partial  volume.  Model  parameters  ( n ,  m, 
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TD50/5)  used  in  this  study  were  those  fitted  to  the  model  by  Burman  etal  (1991)  for  the  normal 
tissue  tolerance  data  compiled  by  Emami  et  al  (1991). 

2.5.  Computational  algorithm 

After  considering  the  automatic  trade-off  strategy,  the  generalized  objective  function  takes  the 
following  form: 

1,  .  n, 

F  =  ^>r—  T  {1  +v[[Dc(i)/Dr,rc(Y'n'}\Dc(i)  -  Z50r(i)|*' 

r=l  T  i=I 

s„  1  N„ 

+  +  <lDc(i)/Drl,rcfY/n'’}Dc(i)k°.  (6) 

<7=1  1  =  1 

We  implemented  a  software  module  to  optimize  the  objective  function  (6)  in  the  PLUNC 
treatment  planning  system  (University  of  North  Carolina,  Chapel  Hill,  NC).  The  ray-by-ray 
iterative  algorithm  (SIITP)  reported  earlier  (Xing  et  al  1998)  was  employed  to  obtain  the 
optimal  beam  intensity  profiles.  The  values  of  kT  and  k„  in  equation  (6)  were  set  to  be  2,  but 
the  behaviour  of  the  system  for  a  few  other  combinations  of  kT  and  k„  were  also  checked  for  the 
prostate  case.  The  reference  dose,  D„  rcf,  was  chosen  to  be  TD5/5  of  the  corresponding  critical 
organ.  For  the  target,  DT  rc f  was  set  as  the  prescription  dose.  Figure  1  shows  the  flow  chart 
of  the  calculation  process.  In  the  current  study  we  specify  a  maximum  number  of  iterations 
as  the  termination  condition  of  the  optimization  process.  The  DVHs  can  be  inspected  in  each 
iterative  step  to  visually  monitor  the  optimization  process. 

2.6.  Case  studies 

Two  cases,  a  prostate  case  and  a  head-neck  case,  were  used  to  evaluate  the  proposed  inverse 
planning  formalism.  The  optimization  results  were  compared  with  those  obtained  using  the 
conventional  dose-based  optimization  method,  which  was  described  in  detail  by  Xing  et  al 
(1998).  The  optimization  parameters  in  the  dose-based  method  were  adjusted  by  trial-and- 
error  to  obtain  an  ‘optimal’  plan. 

In  the  prostate  case,  the  target  volume  included  the  prostate  and  seminal  vesicles.  The 
sensitive  structures  included  rectum,  bladder  and  femoral  heads.  All  the  IMRT  plans  used 
identical  configuration  of  five  equally  spaced  15  MV  photon  beams  with  gantry  angles  of 
0°,  72°,  144°,  216°  and  288°  (in  IEC  convention).  The  plans  were  normalized  to  deliver  the 
prescription  dose  of  70  Gy  to  99%  of  the  target  volume.  The  parameter  nT  was  chosen  to  be 
—0.2  for  the  target.  The  parameters  used  in  the  NTCP  calculations  of  the  rectum,  bladder  and 
femoral  heads  are  listed  in  table  1.  Table  2  summarizes  the  optimization  parameters  for  both 
the  newly  proposed  and  dose-based  approaches. 

For  the  prostate  case,  we  also  studied  the  influence  of  two  more  combinations  of  ka  and 
kT.  These  included  (ka  =  2,  kz  =  4)  and  (k„  =  0,  kx  =  2).  In  the  latter  case,  we  have  included 
a  higher  order  term  of  the  dose-volume  effect  (the  third  term  in  equation  (3)  with  r]2  =  1) 
to  ensure  that  the  high  effective  volume  voxels  are  penalized  enough  in  the  absence  of  the 
dose-based  factor. 

In  the  head-and-neck  case,  the  organs  at  risk  included  the  eyes,  optic  nerves,  optic  chiasm, 
brainstem,  spinal  cord  and  parotids.  Two  targets  were  the  gross  target  volume  (GTV)  and  the 
clinical  target  volume  (CTV),  which  includes  the  microscopic  disease  region  surrounding  the 
GTV.  The  plan  was  normalized  to  deliver  a  prescription  dose  of  70  Gy  to  at  least  99%  of 
the  GTV  and  62  Gy  to  at  least  95%  of  the  CTV.  Nine  equally  spaced  6  MV  coplanar  beams 
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Figure  1.  A  flow  chart  of  the  proposed  optimization  process. 

Table  1.  The  radiological  parameters  for  various  sensitive  structures  used  in  this  study. 


Sensitive  structures  n 

m 

O50/5  (Gy) 

£>s/5  (Gy) 

Bladder 

0.50 

0.11 

80 

65 

Rectum 

0.12 

0.15 

80 

60 

Femoral  head 

0.25 

0.12 

65 

52 

Eye  lens 

0.3 

0.27 

18 

10 

Optic  nerve 

0.25 

0,14 

65 

50 

Optic  chiasm 

0.25 

0.14 

65 

50 

Spinal  cord 

0.05 

0.175 

66.5 

47 

Brainstem 

0.16 

0.14 

65 

50 

Parotid 

0.70 

0.18 

46 

32 

(0°,  40°,  80°,  120°,  160°,  200°,  240°,  280°  and  320")  were  used  for  this  case.  The  parameter 
nT  was  —0.5  for  both  GTV  and  CTV.  The  parameters  used  for  the  computation  of  the  NTCPs 
of  the  sensitive  structures  are  also  obtained  from  the  same  source  stated  earlier  and  are  listed 
in  table  1.  The  optimization  parameters  for  both  the  two  techniques  are  summarized  in 
table  3. 
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Table  2.  Summary  of  the  optimization  parameters  used  in  the  dose-based  and  proposed  approaches 
for  the  prostate  case. 


Organs 

The  dose-based  approach 

The  proposed 
approach 
Generic 
importance 
factors  (r|) 

Relative 

importance 

factors 

Target  prescription 
and  OAR  tolerance 
doses  (Gy) 

Target 

5.0 

78 

5 

Bladder 

1.2 

48 

2 

Rectum 

1.8 

43 

2 

Femoral  head  (R) 

1.0 

32 

1 

Femoral  head  (L) 

1.0 

32 

1 

Normal  tissue 

0.5 

65 

0.3 

Table  3.  Summary  of  the  optimization  parameters  used  in  the  dose-based  and  proposed  approaches 
for  the  head-and-neck  case. 


Organs 

The  dose-based  approach 

The  proposed 
approach 
Generic 
importance 
factors  (rf) 

Relative 

importance 

factors 

Target  prescription 
and  OAR  tolerance 
doses  (Gy) 

GTV 

3.0 

70 

4.0 

CTV 

4.0 

62 

6.0 

Spinal  cord 

2.0 

30 

3.0 

Brainstem 

1.5 

30 

2.0 

Left  optic  nerve 

1.0 

25 

1.0 

Right  optic  nerve 

1.0 

25 

1.0 

Left  eye 

2.0 

6 

3.0 

Right  eye 

2.0 

6 

3.0 

Left  parotid 

1.2 

25 

1.0 

Right  parotid 

1.2 

25 

1.0 

Optic  chiasm 

1.0 

25 

1.0 

Normal  tissue 

0.5 

40 

0.5 

3.  Results 

3.1.  Prostate  1MRT plans 

Figures  2  and  3  summarize  the  results  of  the  two  IMRT  plans  obtained  using  the  newly 
proposed  and  conventional  techniques.  Figure  2  compares  the  isodose  distributions  in  two 
transverse  slices  and  a  sagittal  slice  for  the  two  plans.  The  DVHs  of  the  target  and  sensitive 
structures  are  plotted  in  figure  3,  in  which  the  solid  and  dashed  lines  represent  the  DVHs 
obtained  using  the  new  and  conventional  approaches,  respectively.  The  calculated  NTCPs  of 
rectum,  bladder  and  femoral  heads  for  both  IMRT  plans  are  listed  in  table  4.  According  to  the 
table,  it  is  seen  that  the  NTCPs  of  the  sensitive  structures  are  improved  significantly.  For  the 
rectum,  for  example,  the  NTCP  is  reduced  from  0.45%  to  0.03%.  Our  results  also  indicate  that 
the  main  compromise  in  a  prostate  IMRT  treatment  seems  to  be  between  the  tumour  coverage 
and  the  rectum  complication  because  the  NTCP  of  the  rectum  is  much  higher  than  that  of  other 
sensitive  structures. 
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Axial  1  Axial  2  Sagittal 


(b) 


Figure  2.  Comparison  of  the  isodose  distributions  of  the  two  prostate  IMRT  plans:  (a)  the 
conventional  dose-based  approach;  (b)  the  newly  proposed  approach. 


Table  4.  Comparison  of  the  NTCP  for  the  two  IMRT  plans  for  the  prostate  case. 


NTCP (%) 

The  dose -based  IMRT  plan 

The  proposed  IMRT  plan 

Bladder 

0.017 

0.00030 

Rectum 

0.45 

0.029 

Femoral  head  (R) 

0.000076 

0.0000038 

Femoral  head  (L) 

0.000032 

0.000015 

The  above  results  demonstrate  that,  for  comparable  target  coverage,  the  new  inverse 
planning  technique  greatly  improves  the  critical  structure  sparing,  especially  the  rectum 
sparing.  By  comparing  the  isodose  distributions  of  the  two  plans  (figure  2),  it  is  seen  that  the 
dose  gradient  at  the  interface  between  the  target  and  the  rectum  is  much  steeper  for  the  IMRT 
plan  obtained  with  the  new  formalism.  Furthermore,  it  is  intriguing  that  the  non-sensitive 
structure  normal  tissue  also  receives  fewer  doses  in  comparison  with  that  of  the  dose-based 
optimization.  Our  results  suggest  that  the  improvement  in  the  critical  structure  sparing  is 
achieved  not  at  the  cost  of  higher  target  dose  inhomogeneity,  which  is  commonly  seen  in 
IMRT  plan  optimization. 

The  resultant  DVHs  when  ka  =  2  and  kT  was  increased  from  2  to  4  in  equation  (6)  are 
plotted  in  figure  3  as  the  dotted  curves.  While  the  target  dose  uniformity  is  improved  when 
the  kx  increases,  the  doses  to  the  rectum  and  bladder  are  worsened.  The  results  make  intuitive 
sense  as  when  the  kr  increases,  more  penalty  is  applied  towards  dosimetric  deviation  from  the 
prescription.  The  DVHs  when  the  objective  function  for  the  target  is  a  hybrid  of  dose-volume- 
and  dose-based  functions  ( kx  =  2)  and  that  for  the  sensitive  structures  are  purely  dose-volume 
based  (k„  =  0)  are  shown  in  figure  3  as  dash-dotted  curves.  In  this  case,  a  high  order  term 
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Figure  3.  Comparison  of  DVHs  of  the  prostate  IMRT  plans  obtained  using  the  proposed  approach 
(solid  curves)  and  the  conventional  dose-based  approach  (dashed  curves).  The  dotted  curves 
represent  the  results  obtained  with  kT  =  4  and  k„  =  2  in  equation  (6).  The  dash-dotted  curves  are 
the  DVHs  with  kx=  2  and  kn  =  0  (a  higher  order  term,  the  third  term  in  equation  (3),  was  included 
during  the  optimization). 


of  the  dose-volume  effect  in  equation  (3)  was  added  to  ensure  that  the  high  effective  volume 
voxels  are  penalized  enough  in  the  absence  of  the  dose-based  factor.  Interestingly,  as  can  be 
seen  from  figure  3,  the  results  so  obtained  were  very  similar  to  that  obtained  with  the  hybrid 
objective  function. 

3.2.  Head-and-neck  IMRT  plans 

Figure  4  shows  the  isodose  distributions  obtained  using  the  two  different  planning  techniques 
in  three  transverse  slices,  one  sagittal  slice  and  one  coronal  slice  for  the  two  plans.  Figure  5 
compares  the  DVHs  of  the  targets  and  sensitive  structures,  in  which  the  solid  and  dashed 
lines  represent  the  DVHs  obtained  using  the  newly  proposed  and  conventional  approaches, 
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(a)  Axial  1  (a)  Axial  2  (a)  Axial  3  (a)  Coronal 


(b)  Axial  1  (b)  Axial  2  (b)  Axial  3  (b)  Coronal 


(a)  Sagittal  (b)  Sagittal 


Figure  4.  Comparison  of  the  isodose  distributions  of  the  two  head-and-neck  IMRT  plans:  (a)  the 
conventional  dose-based  approach;  (b)  the  newly  proposed  approach. 

Table  5.  Comparison  of  the  NTCP  for  the  two  IMRT  plans  for  the  head-and-neck  case. 


NTCP  (%) 

The  dose-based  IMRT  plan 

The  proposed  IMRT  plan 

Spinal  cord 

0.043 

0.0025 

Brainstem 

0.012 

0.0040 

Left  Eye 

0.27 

0.18 

Right  Eye 

0.24 

0.12 

Left  parotid 

0.21 

0.056 

Right  parotid 

0.22 

0.064 

Optic  chiasm 

0.00024 

0.000064 

Left  optic  nerve 

0.000064 

0.0000075 

Right  optic  nerve 

0.000043 

0.000025 

respectively.  The  calculated  NTCPs  of  eyes,  optical  nerves,  optical  chiasm,  brainstem,  spinal 
cord  and  parotids  for  both  plans  are  shown  in  table  5.  As  seen  from  the  isodose  distributions 
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(b)  Dose  (Gy) 

Figure  5.  Comparison  of  the  DVHs  of  the  two  head-and-neck  IMRT  plans  obtained  using  our 
newly  proposed  approach  (solid  curves)  and  the  conventional  dose-based  approach  (dashed  curves). 


(figure  4)  and  DVHs  (figure  5),  with  comparable  GTV  and  CTV  dose  coverage  and  dose 
homogeneity,  the  doses  to  the  sensitive  structures  are  dramatically  reduced.  The  dose  reduction 
is  particularly  pronounced  in  the  spinal  cord,  brainstem,  parotids  and  eyes.  For  the  left  and 
right  parotids,  for  example,  the  fractional  volume  receiving  a  dose  above  25  Gy  is  reduced  from 
35%  to  20%  and  1 5%,  respectively.  Consistent  with  the  enhanced  dosimetric  conformality  and 
similar  to  the  prostate  case,  much  steeper  dose  gradient  occurs  near  the  boundary  of  the  target 
volume.  The  dose  to  the  non-sensitive  structure  normal  tissue  is  also  lower  in  comparison 
with  the  conventional  IMRT  plan.  While  the  NTCPs  of  the  optical  nerves  and  optical  chiasm 
are  small  and  it  is  difficult  to  draw  a  conclusion,  from  table  5,  it  is  quite  clear  that  the  NTCPs 
of  the  eyes,  parotids,  spinal  cord  and  brainstem  are  improved  significantly.  We  emphasize 
once  again  that  the  significant  improvement  in  sensitive  structure  sparing  is  achieved  without 
deteriorating  the  dose  coverage  of  the  GTV  and  CTV. 
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(d)  Dose  (Gy) 


Figure  5.  (Continued.) 


4.  Discussion 

The  currently  available  dose-based  objective  functions  do  not  truly  reflect  the  nonlinear 
relationship  between  the  dose  and  the  response  of  tumours  and  tissues,  and  it  is  highly 
desirable  to  incorporate  clinical  outcome  data  in  the  formulation  of  inverse  planning  to  guide 
the  plan  optimization  process.  While  the  dose  dependence  of  a  clinical  endpoint  may  be 
degenerate  in  the  sense  that  it  may  be  caused  by  a  variety  of  dose  distributions  or  DVHs, 
there  exists  no  mechanism  in  conventional  inverse  planning  to  model  the  phenomenon.  The 
irradiation  of  parotid  glands  mentioned  in  the  introduction  represents  an  example  of  this.  On 
the  other  hand,  the  conventional  objective  function  may  impose  some  unrealistic  degeneracy 
that  is  inconsistent  with  clinical  experience.  For  example,  assume  that  in  a  treatment  plan  the 
prescription  dose  to  the  target  is  70  Gy  and  that  the  tumour  is  divided  into  two  parts  with  the 
same  volume,  one  receiving  a  dose  60  Gy  and  another  80  Gy.  The  penalty  values  of  the  two 
scenarios  would  be  the  same  according  to  the  conventional  quadratic  function.  Obviously,  the 
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two  different  dose  distributions  would  lead  to  different  outcomes.  The  cold  part,  which  would 
greatly  diminish  the  tumour  control,  is  more  detrimental  than  the  hot  spot. 

In  this  paper  we  have  established  a  general  inverse  planning  framework  in  which  the 
penalty  at  a  voxel  depends  not  only  on  the  dose  deviation  from  the  desired  value  but  also 
the  dose-volume  status  of  the  involved  organs.  The  technique  circumvents  the  problems 
mentioned  above  and  makes  it  possible  to  take  advantage  of  the  clinical  outcome  data.  Our 
study  shows  that  the  incorporation  of  existing  clinical  knowledge  can  greatly  facilitate  the 
inverse  planning  process  and  allows  us  to  obtain  better  IMRT  plans.  For  the  same  target  dose 
coverage,  the  critical  structure  sparing  was  substantially  improved  for  both  cases.  Physically, 
we  believe  that  the  superior  performance  of  the  new  formalism  arises  from  the  adequate 
modulation  of  the  voxel  dependent  weighting  induced  by  the  dose-volume  factor  /  (see 
equations  (3)  and  (4)).  In  conventional  dose-based  objective  function,  /  =  1,  and  a  tacit 
assumption  that  all  points  within  a  structure  are  equivalent  has  been  made.  The  use  of  dose- 
volume  factor/ given  by  equation  (3)  enables  us  to  weight  different  voxels  according  to  the 
local  doses.  In  this  way,  we  can  effectively  ‘boost’  those  target  regions  where  the  doses  are  low 
or  penalize  more  those  sensitive  structure  regions  where  the  doses  are  high.  The  dose-volume 
induced  voxel  inequality  is  an  important  feature  of  the  new  inverse  planning  formalism  and  is 
the  main  driving-force  in  improving  the  dose  distributions. 

The  use  of  clinical  knowledge  can  also  facilitate  the  determination  of  the  structure  specific 
importance  factors.  While  the  general  influence  of  the  importance  factors  on  the  solution  is 
known,  the  specific  response  of  the  plan  to  a  variation  in  the  factors  is  not  clear  until  the  dose 
optimization  is  done,  which  necessitates  a  manual  trial-and-error  adjustment  of  the  factors 
to  achieve  an  acceptable  trade-off.  The  underlying  deficiency  of  the  conventional  approach 
is  that  the  importance  factors  are  purely  heuristic  and  lack  physical/clinical  meanings.  In 
this  work  we  proposed  a  new  scheme  for  modelling  the  trade-off  and  develop  an  algorithm 
to  auto-detcrmine  the  factors.  The  importance  of  an  organ  was  written  into  a  product  of  a 
generic  and  a  dose-dependent  factor.  The  latter  was  related  to  the  corresponding  TCP  or 
NTCP.  After  the  beam  optimization,  the  dose-dependent  factors  were  increased  or  decreased 
according  to  the  values  of  TCPs/NTCPs.  This  procedure  is  similar  to  that  reported  by  Xing 
et  al  (1999),  where  a  DVH-based  ‘distance’  was  used  for  the  assessment  of  the  trade-off 
status  after  each  optimization.  In  reality,  other  types  of  plan  evaluation  indices,  such  as 
mean/maximum/minimum  doses,  can  also  be  employed  for  the  purpose.  We  noted  that,  with 
the  use  of  a  new  objective  function,  the  final  solution  becomes  much  less  sensitive  to  choice 
of  {r*}.  This  feature  may  have  practical  implications  in  simplifying  the  inverse  planning 
process. 

We  should  acknowledge  that  the  new  technique,  just  like  any  other  dose  response-based 
technique,  may  be  limited  by  the  scarcity  and  uncertainty  of  biological  data  and  the  limited 
predictive  power  of  the  TCP/NTCP  models.  It  is  important  to  note,  however,  that  the 
TCP/NTCP  is  used  as  a  relative  ranking  in  our  plan  optimization  algorithm  instead  of  a 
clinical  decision-making  tool.  Because  of  the  phenomenological  nature  of  the  modelling,  one 
may  further  modify  the  structure  specific  importance  manually  to  achieve  a  certain  clinical 
goal.  The  proposed  technique  can,  at  least,  provide  us  with  a  good  starting  point  for  the 
fine-tuning.  Our  experience,  along  with  the  results  shown  in  the  above  section,  indicates  that 
the  technique  proposed  in  this  work  is  capable  of  generating  clinically  sensible  plans  and  is 
much  more  efficient  than  the  manual  selection  process. 

Finally,  we  mention  that  the  equivalent  uniform  dose  (EUD)-based  dose  optimization 
(Niemierko  1997,  Wu  et  al  2002,  Thieke  et  al  2003,  Lian  and  Xing  2004)  can  be  cast  into  the 
realm  of  the  above  inverse  planning  framework  and  represents  a  special  case  of  the  general 
formalism  described  in  this  work.  Indeed,  assuming  EUD  =  [jj  D“(i))l/a ,  a  =  1  /n  and 
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setting  kr  and  ka  in  equation  (4)  to  zero,  we  can  rewrite  equation  (4)  into 

u 

F  =  EM1  +  V,(EUD/EUDtrcf)a'  +  ,'(EUD/EUD,rcf)^'  +  •  •  •  ] 

r=l 

sa 

+  J2r"['+  tf(EUD/EUD„rcfr»  +  ^2(EUD/EUDCTrcf)2fl"  +  ••■],  (7) 

<7  =  1 

which  becomes  a  function  of  EUD.  Different  from  the  EUD-based  model,  the  general  hybrid 
objective  function  given  in  equation  (6)  treats  the  dose-volume  effect  at  a  more  fundamental 
voxel  level  with  the  actual  radiation  dose  considered,  which  is  more  flexible  than  the  EUD 
defined  at  a  structure  level.  Because  of  this,  other  clinical /dosimetric  requirements  can  easily 
be  integrated.  Our  study  for  the  prostate  case  suggests  that  it  is  necessary  to  include  the  higher 
order  contribution(s)  if  equation  (3)  or  (7)  is  used  to  appropriately  model  a  sensitive  structure. 
Alternatively,  a  hybrid  of  dose-volume  and  dose-based  objective  function,  as  given  by 
equation  (6),  can  yield  equally  good  plans.  In  practice,  equation  (6)  is  quite  broad 
and  seems  to  model  the  inverse  planning  system  effectively.  It  may  also  find  a  natural 
application  in  functional  image-guided  IMRT,  where  the  goal  is  generally  to  produce  a  spatially 
inhomogeneous  dose  distribution  (Xing  et  al  2002).  Finally,  we  note  that  the  formalism  does 
not  involve  the  prescription  of  EUD,  which  could  be  problematic  for  practical  implementation 
of  an  EUD-based  model. 

5.  Conclusion 

Inverse  planning  is  an  important  step  in  IMRT  and  its  performance  crucially  determines  the 
quality  of  IMRT  treatment  plans.  In  this  work,  we  provide  a  mechanism  for  incorporating 
clinical  endpoint  data  into  the  inverse  treatment  planning  process  and  established  a  clinically 
practicable  inverse  planning  framework.  We  employed  the  effective  volume  in  voxel  domain 
to  take  the  dose-volume  effect  of  the  involved  organs  into  account.  The  new  formalism 
gives  important  insight  into  the  problem  of  therapeutic  plan  optimization.  An  algorithm 
for  using  computers  to  aid  the  determination  of  structure  specific  importance  factors  was 
also  developed.  A  key  step  for  accomplishing  the  auto-determination  of  the  importance 
factors  is  the  decomposition  of  the  conventional  importance  factor  into  a  generic  importance 
and  a  dose-dependent  component.  Two  case  studies  were  presented  to  demonstrate  the 
advantages  of  the  proposed  objective  function.  Comparison  of  the  newly  proposed  approach 
with  the  conventional  inverse  planning  technique  indicated  that  the  algorithm  is  capable  of 
greatly  improving  the  sensitive  structure  sparing  with  comparable  target  dose  coverage  and 
homogeneity. 
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Radiobiological  treatment  planning  depends  not  only  on  the  accuracy  of  the  models  describing  the 
dose-response  relation  of  different  tumors  and  normal  tissues  but  also  on  the  accuracy  of  tissue 
specific  radiobiological  parameters  in  these  models.  Whereas  the  general  formalism  remains  the 
same,  different  sets  of  model  parameters  lead  to  different  solutions  and  thus  critically  determine  the 
final  plan.  Here  we  describe  an  inverse  planning  formalism  with  inclusion  of  model  parameter 
uncertainties.  This  is  made  possible  by  using  a  statistical  analysis-based  frameset  developed  by  our 
group.  In  this  formalism,  the  uncertainties  of  model  parameters,  such  as  the  parameter  a  that 
describes  tissue-specific  effect  in  the  equivalent  uniform  dose  (EUD)  model,  are  expressed  by 
probability  density  function  and  are  included  in  the  dose  optimization  process.  We  found  that  the 
final  solution  strongly  depends  on  distribution  functions  of  the  model  parameters.  Considering  that 
currently  available  models  for  computing  biological  effects  of  radiation  are  simplistic,  and  the 
clinical  data  used  to  derive  the  models  are  sparse  and  of  questionable  quality,  the  proposed  tech¬ 
nique  provides  us  with  an  effective  tool  to  minimize  the  effect  caused  by  the  uncertainties  in  a 
statistical  sense.  With  the  incorporation  of  the  uncertainties,  the  technique  has  potential  for  us  to 
maximally  utilize  the  available  radiobiology  knowledge  for  better  IMRT  treatment.  ©  2004 
American  Association  of  Physicists  in  Medicine.  [DOI:  10.1118/1.1785451] 

Key  words:  inverse  planning,  dose  optimization,  biological  models,  IMRT 


INTRODUCTION 

Most  intensity-modulated  radiotherapy  (IMRT)  optimization 
systems  at  present  use  dose  and/or  dose  volume-based  objec¬ 
tive  functions,1-6  which  guide  the  IMRT  planning  by  impos¬ 
ing  a  penalty  according  to  the  difference  between  the  com¬ 
puted  and  prescribed  doses.  A  well-known  drawback  of  the 
dose-based  inverse  planning  is  that  the  nonlinear  dose  re¬ 
sponse  of  tumor  or  normal  structures  is  not  fully  considered. 
A  number  of  mathematical  models  have  been  developed  over 
the  years  to  better  describe  the  biological  effect  of  radiation, 
which  include  tumor  control  probability  (TCP),7  normal  tis¬ 
sue  complication  probability  (NTCP),8  equivalent  uniform 
dose  (EUD)9  and  the  probability  of  uncomplicated  tumor 
control  (/>+). 10,11  In  parallel  to  these  modeling  efforts,  con¬ 
siderable  works  have  also  been  done  to  use  these  biological 
models  to  construct  more  meaningful  objective  functions  for 
therapeutic  dose  optimization.1-16 

Generally  speaking,  radiobiological  formalism  involves 
the  use  of  model  parameters  that  are  of  considerable 
uncertainty.7,17-22  For  instance,  the  radiosensitivity  a  of 
Webb’s  TCP  model  varies  from  0.157  to  0.090  Gy-1  when 
model  parameters  were  fit  to  103  patients’  data.7  Biological 
“margins”  have  been  used  to  account  for  the  variability  in 
radiation  sensitivity.  Similar  to  the  use  of  a  safety  margin  to 
account  for  the  potential  uncertainties  in  targeting  a  tumor, 
this  method  assigns  more  conservative  radiosensitivity  val¬ 
ues  to  the  tumor  or  sensitive  structures  to  deal  with  the  po¬ 
tential  uncertainty  of  the  parameter.23  Kaver  et  al.  proposed  a 
stochastic  optimization  to  account  for  clinical  uncertainties, 


including  the  varying  radiosensitivity.24,25  The  objective 
function  was  constructed  based  on  a  linear  quadratic  Poisson 
model  which  approximates  the  probability  of  curing  the  pa¬ 
tient  or  inflicting  injury.  Two  parameters  in  the  model  could 
be  calculated  if  the  standard  deviation  of  dose  per  faction 
was  known.  The  optimization  was  thus  executed  correspond¬ 
ing  to  different  standard  deviations. 

We  have  recently  presented  a  general  statistical  analysis- 
based  inverse  planning  framework26,27  and  applied  it  to  in¬ 
vestigate  the  influence  of  model  parameter  uncertainties  in 
biologically  based  dose  optimization.28  The  purpose  of  this 
paper  is  to  provide  a  detailed  description  of  the  technique 
and  addresses  several  important  issues  related  to  the  dose 
optimization  in  the  presence  of  model  parameter  uncertain¬ 
ties.  In  our  approach,  the  uncertainty  of  a  model  parameter  is 
quantified  by  a  probability  density  function  and  its  influence 
is  then  incorporated  into  inverse  planning  through  the  use  of 
a  statistical  inference  theorem.26  The  technique  is  illustrated 
by  using  a  hypothetical  C-shaped  tumor  case,  a  prostate  tu¬ 
mor  case  and  a  paraspinal  tumor  case  with  an  EUD-based 
model.  Considering  that  currently  available  models  for  com¬ 
puting  biological  effects  of  radiation  are  simplistic,  and  the 
data  they  rely  on  are  sparse  and  of  questionable  quality,  the 
proposed  technique  provides  us  with  an  effective  tool  to 
minimize  the  effect  caused  by  the  uncertainties  in  a  statistical 
sense.  The  treatment  plans  so  obtained  are  generally  less 
sensitive  to  the  inter-patient  variation  and  other  types  of  un¬ 
certainties  that  may  otherwise  influence  the  final  treatment 
plan  greatly. 
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METHODS  AND  MATERIALS 
Statistical  analysis-based  inverse  planning 

The  inverse  problem  as  posed  for  IMRT  consists  of  the 
determination  of  the  beamlet  weight  vector  w  when  a  desired 
plan  is  prescribed.  In  a  vectorial  form,  the  dose  to  the  points 
in  the  treatment  region  depends  upon  the  beamlet  weights  w 
as 

Dc  —  d'\v,  (1) 

where  d  represents  the  dose  deposition  coefficients  matrix, 
expressing  the  dose  deposited  to  any  patient  point  when  ir¬ 
radiated  with  a  unit  weight  beamlet.  The  total  number  of 
physically  realizable  dose  distributions  Dc  in  IMRT  is  enor¬ 
mous  and  increases  exponentially  with  the  number  of  beam- 
lets.  Inverse  planning  is  essentially  a  plan  selection  process 
from  the  vast  pool  of  physically  realizable  solutions.  In  a 
recent  paper,  Xing  et  al.26  introduced  a  statistical  analysis- 
based  inverse  planning  technique.  In  this  approach  the  com¬ 
monly  used  objective  function  is  reformulated  into  a  prob¬ 
ability  density  function  whose  value  gives  the  figure  of  merit 
of  a  dose  distribution.  A  virtue  of  the  approach  is  that  it 
allows  us  to  obtain  solution  in  the  presence  of  uncertainties 
of  the  prescription  parameters  or  other  model  parameters  us¬ 
ing  a  statistical  inference  technique.  Application  of  the  tech¬ 
nique  to  deal  with  a  system  with  a  set  of  variable  dose  pre¬ 
scriptions  has  been  described  in  another  work  of  our  group.27 
Here  we  use  the  formalism  for  biological  modeling  based- 
inverse  planning  in  the  presence  of  model  parameter  uncer¬ 
tainties.  To  be  specific,  we  use  an  equivalent  uniform  dose 
(EUD)-based  objective  function  employed  by  Wu  et  al. 13,29 
and  discuss  the  consequences  of  the  variation  of  model  pa¬ 
rameter  a  and  how  to  incorporate  the  fluctuations  into  in¬ 
verse  planning  dose  optimization  to  obtain  statistically  opti¬ 
mal  solutions. 


Fig.  1.  The  flow  chart  of  the  optimization  process  with  the  inclusion  of 
model  parameter  uncertainty. 


0’ 


10  2D  30  40  50  BO 


Fig.  2.  (A)  The  sketch  of  the  hypothetical  case  with  C-shaped  target  and  the 
beam  setup  for  dose  optimization.  (B)  The  dose  distribution  corresponding 
to  the  parameters  listed  in  Table  1  and  the  probabilistic  distribution  shown  in 
Fig.  3  d2. 

EUD  model  and  EUD-based  objective  function 

The  concept  of  equivalent  uniform  dose  (EUD)  for  tumor 
was  originally  introduced  by  Niemierko  as  the  biologically 
equivalent  dose  that,  if  given  uniformly,  would  lead  to  the 
same  cell  kill  in  the  tumor  volume  as  the  actual  nonumiform 
dose  distribution.  Recently,  Niemierko  et  al.  suggested  a 
phenomenological  form9, 13,30 
/  j  \  i/o 

EUD=(-E  D"J  ,  (2) 

for  both  tumor  and  normal  tissues,  where  N  is  the  number  of 
voxels  in  the  structure,  Dt  is  the  dose  delivered  to  the  /th 
voxel,  a  is  the  tumor  or  normal  tissue-specific  parameter  that 
describes  the  dose-volume  effect.  EUD  described  in  Eq.  (2) 
is  the  general  mean  of  the  non-uniform  dose  distribution. 
According  to  the  mathematic  properties  of  the  function,31 
for  o  =  oo  j  the  EUD  is  equal  to  the  maximum  dose,  and  for 
a  =—oo,  the  EUD  is  equal  to  the  minimum  dose.  Tumors 
generally  have  large  negative  values  of  a,  whereas  serial 
critical  structures  (e.g.,  spinal  cord  and  rectum)  have  large 
positive  values  and  parallel  critical  structures  that  exhibit  a 
large  dose-volume  effect  (e.g.,  liver,  parotids,  and  lungs) 
have  small  positive  values. 
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Table  I.  The  conventional  EUD-based  optimization  parameter  for  the  hy-  Table  II.  The  conventional  EUD-based  optimization  parameter  for  prostate 
pothetical  IMRT  treatment  of  a  C-shaped  tumor,  cancer. 


"Contains  parameters  for  the  target  treated  as  virtual  normal  tissue  to  limit  "Contains  parameters  for  the  target  treated  as  virtual  normal  tissue  to  limit 
dose  inhomogeniety.  dose  inhomogeniety. 


The  objective  function  or  figure  of  merit  used  to  measure 
the  goodness  of  a  dose  distribution  and  guide  the 
optimization.1  In  the  present  paper,  the  system  objective 
function  is  given  by13 


F=Ufj, 

j 

where  the  component  subcore  fj  may  be  either 

1 

fr~  /EUD0\”’ 

\  eud) 

for  tumors,  or 


(3) 


(4) 


(5) 


for  normal  tissues  and  organs  at  risk  (OARs).  EUD0  is  the 
desired  dose  parameter  for  the  target  volume  and  the  maxi¬ 
mum  tolerable  uniform  dose  for  normal  structures.  Parameter 
n  is  akin  to  the  structure  specific  importance  factor32  in  the 
conventional  inverse  planning  formalism  that  parameterizes 
our  tradeoff  strategy  of  different  structure.  The  large  n  indi¬ 
cates  high  importance. 


Fig.  3.  The  target  and  OAR  DVHs  of  four  optimal  plans  when  parameter  a  is  a  fixed  value  (bar  chart  dl)  and  varies  according  to  three  different  probabilistic 
distributions  (bar  chart  d2,  d3,  and  d4). 
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Fig.  4.  The  EUD  of  the  target  and  objective  function  when  parameter  a  is 
prescribed  according  to  Fig.  3. 


Incorporation  of  the  variation  distribution 
of  the  model  parameter  into  inverse  planning 

We  assume  that  ak  in  the  EUD  model  varies  according  to 
a  simple  Gaussian  distribution 

Pn(ak)  =  P'on  exP{ ~r„[ak~ a0]2},  (6) 

where  a0  is  the  mean  value,  P'0„  is  a  normalization  constant 
and  ak  is  one  of  the  sampling  values  of  a.  For  a  given  dis¬ 
tribution,  the  EUD  and  the  corresponding  figure  of  merit  of 
an  IMRT  plan  vary  with  the  sampling  of  a.  We  thus  rewrite 
Eqs.  (4)  and  (5)  as  conditional  probabilities  for  a  sampled 
ak: 

/V(EUD|a*)= - 

1  + 


EUD0\”’ 
EUD  ) 


P  OAR(EUDk)- 


1 


1  + 


I  EUD 
EUDn 


(8) 


The  objective  function  for  a  structure  m  in  the  presence  of 
uncertainty  in  a  is  expressed  as  the  summation  of  a  series  of 
joint  probabilities 


P„,(EUD)  =  2  Pm(EUD\ak)-PJak),  (9) 

k 


and  the  overall  objective  function  P  of  the  system  is  a  prod¬ 
uct  of  Pm(EUD)  defined  in  Eq.  (9).  That  is 

F=  ln(  1/P)  =  —  In n  -Pm(EUD) 


=  “2  In 2  Pm{VVV\ak)-Pm(ak).  (10) 

m  k 


Optimization  method 

As  described  above,  the  uncertainties  of  model  param¬ 
eters,  {«;},  are  described  by  probability  density  functions 
and  they  are  incorporated  into  the  overall  objective  function 
of  the  system  through  the  joint  probability  given  by  Eq.  (9). 
To  obtain  the  optimal  solution  in  the  presence  of  model  pa¬ 
rameters,  all  we  need  to  do  is  to  minimize  the  overall  objec¬ 
tive  function  given  by  Eq.  (10). 

The  calculation  process  is  schematically  shown  in  Fig.  1 . 
For  the  computational  purpose,  the  probability  density  func¬ 
tion  for  each  structure  is  discreted  into  seven  equally  spaced 
points.  We  use  the  Fletcher-Reeves  conjugate  gradient  opti¬ 
mization  algorithm33  to  optimize  the  system.  But  any  other 
iterative  or  stochastic  optimization  can  be  also  employed  to 
optimize  the  system.  A  common  step  in  all  optimization  al¬ 
gorithms  is  the  evaluation  of  the  objective  function  for  a  trial 
beam  profiles  (or  computed  dose  distribution),  which  is 
somewhat  tedious  here  because  of  the  appearance  of  multiple 
ak  s  of  the  involved  structures.  Briefly,  for  a  given  trial  beam 
profiles  or  dose  distribution,  the  evaluation  of  the  objective 
function  consists  of  four  steps:  (i)  For  a  structure  m,  calculate 
the  EUD  corresponding  to  each  possible  ak ;  (ii)  calculate  the 
conditional  probability  for  the  target  and  OAR  using  Eqs.  (7) 
and  (8),  respectively;  (iii)  sum  over  all  possible  ak  to  obtain 
the  joint  probability,  given  by  Eq.  (9);  and  (iv)  sum  over  all 
structures  to  obtain  the  overall  objective  function  value.  Af¬ 
ter  the  dose  optimization,  a  set  of  optimal  beam  profiles  and 
the  corresponding  dose  distribution  and  other  plan  indices 
are  provided  for  the  planner  to  assess  the  clinical  relevance 
of  the  obtained  treatment  plan. 


Test  cases 

The  new  algorithm  was  tested  using  a  hypothetical  phan¬ 
tom  case  with  a  C-shaped  target  and  two  clinical  cases  (a 
prostate  case  and  a  paraspinal  tumor  boost  treatment).  The 
size  of  the  pencil  beam  defined  at  the  isocenter  was  0.5  cm. 
The  configuration  of  the  C-shaped  tumor  case  is  shown  in 
Fig.  2(A).  Nine  6  MV  equi-spaced  beams  were  used  for  the 
treatment  (0°,  40°,  80°,  120°,  160°,  200°,  240°,  280°,  and 
320° — respecting  the  International  Electrotechnical  Commis¬ 
sion  (IEC)  convention).  The  values  of  n  and  a  in  the  EUD- 
based  objective  function  are  listed  in  Table  I.  The  parameter 
a  in  EUD  model  characterizes  the  dose-volume  effect  but  its 
value  is  generally  not  known  accurately  even  for  clinically 
well  studied  organs.  The  influence  of  the  uncertainty  in  the  a 
value  of  a  target  or  sensitive  structure  to  the  final  treatment 
plan  was  studied  and  analyzed. 
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OAR  parameter  distribution 


Fig.  5.  The  target  and  OAR  DVHs  of  four  optimal  plans  when  parameter  a  is  a  fixed  value  (bar  chart  dl)  and  varies  according  to  three  different  probabilistic 
distributions  (bar  chart  d2,  d3,  and  d4). 


Similar  study  was  carried  out  for  the  two  clinical  cases. 
The  six  6  MV  beam  angles  used  for  the  IMRT  prostate  treat¬ 
ment  were  0°,  55°,  135°,  180°,  225°,  and  305°.  Table  II  lists 
some  relevant  parameters  used  for  planning  the  case.  For  the 
IMRT  paraspinal  boost  treatment,  five  6  MV  nonequally 
spaced  coplanar  beams  were  placed  at  the  following  angular 
positions:  95°,  140°,  175°,  225°,  and  275°.  The  target  boost 
dose  was  prescribed  to  16  Gy.  Relevant  parameters  are  listed 
in  Table  II.  The  planning  goal  was  to  find  a  dose  distribution 
that  covered  the  tumor  volume  as  uniformly  as  possible, 
while  maximally  sparing  the  spinal  cord,  liver,  and  kidney. 

RESULTS  AND  DISCUSSIONS 
The  C-shaped  tumor  case 

We  first  investigated  the  behavior  of  the  system  when  the 
parameter  a  of  the  target  EUD  takes  four  different  distribu¬ 
tions,  as  depicted  in  the  bar  charts  shown  on  the  right  of  Fig. 
3,  while  keeping  the  parameter  a  of  the  OAR  at  a  constant 
a0  =  6.0.  In  the  case  shown  in  Fig.  3  dl,  the  parameter  a 
takes  only  a  single  value,  a0=  —  10,  which  is  a  simple  case 
studied  by  Wu  et  a/.13  The  optimal  plans  for  the  four  distri¬ 
butions  of  parameter  a  differ  significantly,  as  indicated  by  the 
target  and  OAR  DVHs  shown  in  Figs.  3(A)  and  3(B).  The 
isodose  plot  corresponding  to  the  a-distribution  shown  in 
Fig.  3  d2  is  plotted  in  Fig.  2(B). 

To  estimate  the  degree  of  sensitivity  of  the  solutions 
against  a  variation  in  a,  we  computed  the  target  EUD  and  the 
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objective  function,  fT,  as  a  function  of  parameter  a  for  the 
four  optimal  dose  distributions  under  different  types  of  un¬ 
certainty  distributions.  The  results  are  plotted  in  Fig.  4.  For 
plan  dl,  the  EUD  changes  from  65  to  71  Gy  when  a  is  varied 
from  - 10  to  —70  and  to  79  Gy  when  a  is  equal  to  140.  The 
objective  function  varies  from  0.11  to  0.85  in  the  range  of 
variation  in  a.  For  plan  d4,  the  EUD  is  narrowed  to  a  range 
between  70  and  79  Gy.  The  EUD  variations  of  plans  d2  and 
d3  are  similarly  reduced.  These  results  suggest  that  the  EUD 
becomes  much  less  sensitive  to  the  variation  in  parameter  a 
in  the  plans  obtained  with  some  “built-in”  distributions  in 
parameter  a  (i.e.,  plans  corresponding  to  Figs.  3  d2  to  d4). 

The  uncertainty  of  parameter  a  of  the  OAR  can  be  simi¬ 
larly  included  in  the  dose  optimization  process  when  its  dis¬ 
tribution  is  known.  In  the  second  study,  we  fixed  the  target 
EUD  parameter  a=  —  10  and  allowed  the  parameter  a  of  the 
OAR  to  take  four  different  distributions  as  plotted  in  the 
right  of  Fig.  5.  The  target  and  OAR  DVHs  for  the  four  pos¬ 
sible  scenarios  are  shown  in  A  and  B.  Once  again,  we  found 
that  the  final  solution  strongly  depends  on  the  distributions  of 
the  parameter  a. 

The  maximum  doses  of  the  OAR  of  the  four  plans  vary 
from  24  to  30  Gy.  Note  that  the  doses  to  the  OAR  in  plans 
d2,  d3,  and  d4  are  less  than  that  of  plan  A,  where  the  param¬ 
eter  a  is  restricted  to  a  single  value,  a0  =  6.  This  is  explain¬ 
able  since  the  parameters  a  in  plans  d2,  d3,  and  d4  are  shifted 
up  to  higher  values.  As  a  increases,  the  EUD  puts  more  em- 
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Fig.  6.  The  EUD  of  the  OAR  and  objective  function  when  parameter  a  is 
prescribed  according  to  Fig.  5. 


phasis  on  the  high  dose  (recall  that  EUD  becomes  the  maxi¬ 
mum  dose  when  a  =  °°).  As  a  consequence  of  the  increased 
“effective”  a  value  in  the  distributions  shown  in  Fig.  5,  d2, 
d3,  and  d4,  the  OAR  dose  is  improved  in  comparison  with 
the  plan  obtained  under  the  assumption  of  a  fixed  a  value 
(Fig.  5,  dl).  Interestingly,  the  target  DVHs  shows  that  four 
distinct  plans  have  very  similar  target  coverage.  It  is  well 
known  that  in  dose  optimization  there  is  generally  no  net 
gain:  an  improvement  in  the  dose  to  a  structure  is  often  ac¬ 
companied  by  a  dosimetrically  adverse  effect(s)  at  other 
points  in  the  same  or  different  structures.  The  result  here 
suggests  that,  from  a  clinical  point  of  view,  it  is  possible  to 
have  a  great  gain  in  one  structure  with  a  little  sacrifice  in 
another  structure.  How  to  find  the  truly  optimal  tradeoff  rep¬ 
resents  a  practical  subject  that  is  worth  of  studying  in  the 
future. 

As  can  be  expected  from  the  discussion  in  previous  para¬ 
graphs,  the  solution  obtained  with  a0  =  6  (Fig.  5,  dl)  is  more 
sensitive  to  a  variation  in  parameter  a.  Indeed,  as  seen  from 
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Fig.  7.  A  transverse  slice  showing  the  anatomical  structures  delineated  for 
the  prostate  tumor  (A)  and  the  corresponding  optimized  dose  distribution 
with  the  parameters  listed  in  Table  11  and  the  probabilistic  distribution 
shown  in  Fig.  8(B). 

Fig.  6,  the  EUD  for  this  plan  varies  from  1  to  30  Gy  when  a 
is  changed  from  —80  to  140.  On  the  other  hand,  the  EUD 
changes  for  the  rest  three  situations  are  much  less  for  the 
same  variation  in  a.  The  upper  bound  of  the  EUD  is  reduced 
to  26  Gy  for  plan  d4,  24  Gy  for  plan  d2,  and  23  Gy  for  plan 
d3.  The  objective  functions  of  four  plans  show  a  similar 
trend. 

The  prostate  tumor  case 

Four  IMRT  plans  with  different  types  of  pre-assumed  un¬ 
certainties  were  generated  for  a  prostate  tumor  case  [Fig. 
7(A)],  These  include:  (i)  The  ^-parameters  for  both  prostate 
target  and  OARs  are  restricted  to  single  values  as  listed  in 
Table  II.  This  plan  serves  as  a  reference  whose  DVHs  are 
shown  in  Figs.  8(A)-8(C)  as  dotted  curves;  (ii)  Only  the 
a-parameter  of  the  prostate  target  takes  a  range  of  values,  as 
depicted  in  the  right  of  Fig.  8(A);  (iii)  Only  the  a-parameter 
of  the  rectum  takes  a  range  of  values,  as  depicted  in  the  right 
of  Fig.  8(B);  and  (iv)  The  a-parameters  of  both  prostate  tar¬ 
get  and  the  rectum  were  allowed  to  take  a  range  of  values,  as 
depicted  in  the  right  of  Fig.  8(C). 
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Fig.  8.  DVHs  for  a  prostate  cancer  case  using  the  conventional  optimization  with  fixed  o-value  (dotted  line)  and  the  newly  proposed  approach  with  the 
inclusion  of  model  parameter  uncertainty  (solid  line).  (A)  Only  the  a-parameter  for  the  target  is  assigned  with  a  probabilistic  distribution;  (B)  Only  the 
a-parameter  for  the  OAR  is  assigned  with  a  probabilistic  distribution;  (C)  Uncertainties  in  the  a-parameter  are  introduced  for  both  the  target  and  OAR. 


DVHs  for  the  plan  using  parameters  defined  in  Table  II 
are  plotted  with  dotted  curves  and  plans  with  the  inclusion  of 
parameter  uncertainty  are  drawn  with  solid  curves  (Fig.  8). 


When  the  parameter  a  in  target  EUD  takes  a  Poisson  distri¬ 
bution  as  shown  in  the  bar  chart  of  Fig.  8(A),  prostate  dose 
homogeneity  is  significantly  improved.  The  minimum  dose 
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Fig.  9.  (A)  A  transverse  slice  showing  the  anatomical  structures  for  a  paraspinal  case;  (B)  DVHs  when  the  a-parameter  for  the  target  is  assigned  with  a 
probabilistic  distribution;  (C)  DVHs  when  the  a-parameter  for  the  OAR  is  assigned  with  a  probabilistic  distribution. 


increases  from  55  to  67  Gy,  and  the  maxim  dose  decreases 
slightly  from  82  to  80  Gy.  However  the  volumes  receiving 
radiation  dose  for  rectum,  bladder  and  normal  tissue  all  in¬ 


crease  significantly  though  the  maximum  dose  remains  simi¬ 
lar.  The  improvement  of  the  target  coverage  and  compromise 
of  OAR  sparing  is  a  natural  outcome  of  the  competitive  re- 
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Table  III.  The  conventional  EUD-based  optimization  parameter  for  paraspi- 
nal  tumor. 


PTV 

PTV" 

Spine 

Liver 

Kidney 

a 

-10.0 

10.0 

6.0 

6.0 

6.0 

EUD0  (Gy) 

16 

17 

12 

6.4 

4.8 

n 

20 

20 

6 

6 

6 

“Contains  parameters  for  the  target  treated  as  virtual  normal  tissue  to  limit 
dose  inhomogeniety. 


quirements  for  targets  and  OARs  imposed  on  the  system. 
The  corresponding  dose  distribution  with  the  target  param¬ 
eter  defined  in  the  bar  chart  A  is  shown  in  Fig.  7(B). 

Next  we  considered  the  inclusion  of  parameter  a  uncer¬ 
tainty  in  EUD  calculation  in  one  of  the  critical  structures- 
rectum  [Fig.  8(B)],  The  irradiated  rectum  volume  for  a  dose 
below  60  Gy  is  less  than  that  of  a  conventional  plan  with  the 
parameter  a  fixed  at  24.  DVHs  for  the  bladder,  normal  tissue 
and  prostate  do  not  change  significantly  compared  to  the  plan 
without  inclusion  of  parameter  uncertainty. 

Lastly,  we  simultaneously  replaced  target  and  rectum  pa¬ 
rameters  with  the  distributions  shown  in  Fig.  8(C).  Similar  to 
that  corresponds  the  prescription  of  Fig.  8(A),  the  prostate 
coverage  is  improved.  However,  the  rectum  DVH  in  this  case 
is  not  worsen  greatly  because  parameter  a  of  rectum  EUD 
was  allowed  to  take  a  spectrum  of  values.  For  bladder  and 
normal  tissue,  although  their  irradiated  volumes  in  the  low 
dose  region  are  higher  than  those  of  the  conventional  plan, 
the  volumes  receiving  high  doses  are  reduced. 

The  paraspinal  tumor  case 

Three  IMRT  plans  were  generated  for  a  paraspinal  tumor 
case  [Fig.  9(A)].  These  include:  (i)  The  a-parameters  for  both 
target  and  OARs  are  restricted  to  single  values  as  listed  in 
Table  III.  This  plan  serves  as  a  reference  whose  DVHs  are 
shown  in  Figs.  9(B)  and  9(C)  as  dotted  curves;  (ii)  Only  the 
a-parameter  of  the  target  takes  a  range  of  values,  as  depicted 
in  the  right  of  Fig.  9(B);  and  (iii)  Only  the  a-parameter  of  the 
spinal  cord  takes  a  range  of  values,  as  depicted  in  the  right  of 
Fig.  9(C). 

When  a  in  target  EUD  takes  a  Poisson  distribution  as 
shown  in  the  bar  chart  of  Fig.  9(B),  dose  homogeneity  is 
slightly  improved.  However,  this  is  achieved  at  the  expense 
of  more  irradiation  to  the  cord.  The  inclusion  of  parameter  a 
uncertainty  in  EUD  calculation  in  the  spinal  cord  [Fig.  9(C)] 
reduced  the  maximum  cord  dose  by  100  cGy.  The  DVHs  of 
the  target,  kidney,  and  liver  were  not  changed  significantly 
compared  to  the  plan  without  inclusion  of  parameter  uncer¬ 
tainty. 

The  influence  of  various  uncertainties  on  the  patient  treat¬ 
ment  has  been  a  subject  of  intense  study.  Fenwick  and 
Nahum  have  included  the  model  parameter  uncertainty  with 
a  standard  deviation  when  calculating  the  NTCP  of 
rectum.34,35  Similarly,  the  inclusion  of  uncertainties  in  the 
patient  setup  and  dose  calculation  has  also  been 
demonstrated.36-38  Deasy  et  al.  have  used  a  bootstrap-based 
method  to  estimate  the  influence  of  biological  parameter  un¬ 


certainties  on  predicting  long-term  salivary  function.21  The 
statistical  method  proposed  here  provides  a  general  frame¬ 
work  to  include  various  uncertainties  in  the  dose  optimiza¬ 
tion  process.  With  minor  modification,  the  technique  can  be 
extended  to  derive  statistically  optimal  solutions  in  the  pres¬ 
ence  of  other  types  of  uncertainties. 

As  can  be  intuitively  imagined,  the  inclusion  of 
a-distribution  will  definitely  change  the  final  dose.  Whether 
it  will  improve  or  worsen  the  final  dose  distribution  will 
generally  depend  on  the  specific  form  of  the  a-distribution, 
and  also  the  metric  used  to  judge  the  goodness  of  a  plan.  If 
the  original  EUD-based  objective  function  is  used  as  the  sole 
metric  for  the  judgment,  the  inclusion  of  a-distribution  may 
make  the  plan  worse.  However,  clinical  decision-making  is 
not  made  by  a  single  function  and  a  “worse”  plan  judged  by 
the  EUD-objective  function  may  turn  out  to  be  clinically 
more  favorable.  In  other  words,  there  is  a  gap  between  math¬ 
ematical  dose  optimization  and  clinical  decision-making. 
The  study  seems  to  suggest  that,  while  it  is  generally  true 
that  there  is  no  net  gain  in  dose  optimization,27  it  is  important 
to  develop  a  method  that  is  capable  of  optimizing  not  only 
the  objective  function  but  also  the  next  level  of  decision¬ 
making.  This  kind  of  optimization  will  allow  us  to  find  the 
solution  that  may  sacrifice  a  little  (i.e.,  clinically  insignifi¬ 
cant)  in  one  or  a  few  structures  but  gain  a  lot  in  other  struc¬ 
tures. 

CONCLUSIONS 

We  have  proposed  and  implemented  a  technique  for  in¬ 
corporating  biological  model  parameter  uncertainties  into  in¬ 
verse  treatment  planning.  The  formalism  is  quite  general  and 
does  not  prerequisite  the  specific  form  of  uncertainty  distri¬ 
butions  of  the  involved  model  parameters.  By  including 
model  parameter  uncertainties,  the  final  solution  becomes 
more  robust  and  the  treatment  outcome  will  be  less  likely 
influenced  by  inter-patient  variation  of  biological  character¬ 
istics.  With  the  increasing  interest  in  radiation  therapy  com¬ 
munity  to  use  biologically  based  models  for  treatment  plan¬ 
ning,  this  work  provides  an  effective  way  to  better  account 
for  the  known  uncertainties  in  the  model  parameters  and  al¬ 
lows  us  to  maximally  utilize  the  available  radiobiology 
knowledge  to  facilitate  patient  care. 
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The  dose  optimization  in  inverse  planning  is  realized  under  the  guidance  of  an  objective  function. 

The  prescription  doses  in  a  conventional  approach  are  usually  rigid  values,  defining  in  most  in¬ 
stances  an  ill-conditioned  optimization  problem.  In  this  work,  we  propose  a  more  general  dose 
optimization  scheme  based  on  a  statistical  formalism  [Xing  et  al.,  Med.  Phys.  21,  2348-2358 
(1999)].  Instead  of  a  rigid  dose,  the  prescription  to  a  structure  is  specified  by  a  preference  function, 
which  describes  the  user’s  preference  over  other  doses  in  case  the  most  desired  dose  is  not  attain¬ 
able.  The  variation  range  of  the  prescription  dose  and  the  shape  of  the  preference  function  are 
predesigned  by  the  user  based  on  prior  clinical  experience.  Consequently,  during  the  iterative 
optimization  process,  the  prescription  dose  is  allowed  to  deviate,  with  a  certain  preference  level, 
from  the  most  desired  dose.  By  not  restricting  the  prescription  dose  to  a  fixed  value,  the  optimiza¬ 
tion  problem  becomes  less  ill-defined.  The  conventional  inverse  planning  algorithm  represents  a 
special  case  of  the  new  formalism.  An  iterative  dose  optimization  algorithm  is  used  to  optimize  the 
system.  The  performance  of  the  proposed  technique  is  systematically  studied  using  a  hypothetical 
C-shaped  tumor  with  an  abutting  circular  critical  structure  and  a  prostate  case.  It  is  shown  that  the 
final  dose  distribution  can  be  manipulated  flexibly  by  tuning  the  shape  of  the  preference  function 
and  that  using  a  preference  function  can  lead  to  optimized  dose  distributions  in  accordance  with  the 
planner’s  specification.  The  proposed  framework  offers  an  effective  mechanism  to  formalize  the 
planner’s  priorities  over  different  possible  clinical  scenarios  and  incorporate  them  into  dose  opti¬ 
mization.  The  enhanced  control  over  the  final  plan  may  greatly  facilitate  the  IMRT  treatment 
planning  process.  ©  2003  American  Association  of  Physicists  in  Medicine. 
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I.  INTRODUCTION 

Inverse  planning  is  used  in  intensity  modulated  radiation 
therapy  (IMRT)  for  deriving  the  optimal  beam  intensity  pro¬ 
files  that  produce  the  best  possible  dose  distribution  for  a 
given  patient.1-16  The  dose  optimization  process  is  usually 
performed  under  the  guidance  of  an  objective  function, 
which  measures  the  “distance”  between  the  physical  and  the 
prescribed  dose  distributions.8,17-20  One  of  the  common  ob¬ 
jective  functions  for  inverse  planning  is  the  quadratic  objec¬ 
tive  function,3,21,22  with  importance  factors  assigned  to  the 
involved  structures  to  prioritize  their  relative  importance  dur¬ 
ing  the  optimization  process.23-25  The  objective  function  is 
defined  as  a  global  quantity  based  on  general  physical  con¬ 
siderations.  When  the  desired  dose  distribution  is  not  attain¬ 
able  during  optimization,  a  compromise  solution  is  found 
using  the  algorithm’s  ranking.  The  compromise  dose  distri¬ 
bution,  however,  is  often  not  what  the  planner  wants  and 
multiple  trial  and  errors  are  needed  to  obtain  a  clinically 
acceptable  IMRT  plan. 

A  main  problem  of  the  existing  IMRT  planning  algo¬ 
rithms  is  the  lack  of  an  effective  mechanism  for  incorporat¬ 
ing  prior  knowledge  into  inverse  planning.31  In  the  past, 
there  have  been  many  attempts  to  introduce  soft/hard  con¬ 
straints  to  steer  the  dose  optimization  process  toward  the 


clinically  desired  solutions.26-30,36  However,  the  constraints 
are  introduced  in  an  ad  hoc  fashion  and  do  not  fully  utilize 
the  partial  information  available  from  years  of  clinical  inves¬ 
tigations  because  of  their  phenomenological  nature.  On  a 
more  fundamental  level,  the  constraints  are  imposed  a  pos¬ 
teriori  and  controls  the  optimization  passively.  Our  purpose 
in  this  paper  is  to  develop  a  statistical  analysis-based  inverse 
planning  formalism  to  more  effectively  utilize  the  prior 
knowledge.  Instead  of  specifying  a  rigid  prescription  dose, 
the  formalism  allows  us  to  use  a  dose  distribution  as  the 
input  prescription  to  the  system,  providing  a  natural  way  for 
us  to  take  advantage  of  the  existing  information  of  the  sys¬ 
tem  variables  and  promising  to  make  the  optimization  out¬ 
come  more  predictable  and  controllable. 

In  the  next  section  we  present  the  details  of  the  new  dose 
optimization  algorithm  after  a  brief  introduction  of  the  con¬ 
cept  of  preference  function.  The  formalism  is  then  applied  to 
a  synthetic  phantom  case  with  C-shaped  tumor  target  and  a 
prostate  case.  Our  results  indicate  that  the  statistical  analysis- 
based  formalism  provides  a  general  framework  for  inverse 
planning  and  is  capable  of  producing  conformal  IMRT  dose 
distribution.  Coupled  with  the  capability  of  the  preference 
function  in  customizing/formalizing  our  prior  clinical  knowl¬ 
edge,  it  is  expected  that  the  proposed  technique  will  have  a 
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broad  implication  and  potential  to  greatly  facilitate  an  IMRT 
planning  process. 

II.  MATERIAL  AND  METHODS 

A.  Theoretical  background 

In  a  vectorial  form,  the  dose  to  the  points  in  the  treatment 
region  depend  upon  the  beamlet  weights  w  as 

Dc=d-w,  (1) 

where  d  represents  the  dose  deposition  matrix,  expressing 
the  dose  deposited  to  any  point  in  the  patient  when  irradiated 
with  a  unit  weight  beamlet  vector.  The  inverse  problem  as 
posed  for  IMRT  is  to  find  a  set  of  beamlet  weights  that  pro¬ 
duce  the  optimal  dose  distribution  by  minimizing  a  therapeu¬ 
tic  objective  function.  The  most  used  objective  function  has  a 
quadratic  form  and  reads32  as 

1  N 

^=77  2  ra[Dc(n)-Dp{n)f,  (2) 

/V  w  =  1 

where  N  is  the  total  number  of  voxels,  r  a  is  the  importance 
factor  that  controls  the  relative  importance  of  a  structure  <r, 
and  Dp  and  Dc  are  prescribed  and  calculated  doses,  respec¬ 
tively. 

In  inverse  planning  algorithm  based  on  the  quadratic  ob¬ 
jective  function  [Eq.  (2)],  the  dose  prescription  to  the  target 
or  sensitive  structure  takes  a  rigid  value.  The  minimization  of 
the  objective  function  is  realized  by  various  algorithms  like 
simulated  annealing,  gradient  methods,  etc. 17,26,32,33  Indepen¬ 
dent  of  the  used  dose  optimization  algorithms,  we  will  call 
these  methods  throughout  the  text  conventional  IMRT  opti¬ 
mization  procedures.  The  problem  is  usually  ill-posed  and 
may  lead  to  negative  fluence  unless  hard  constraints  are 
introduced.3  Practically,  it  is  not  uncommon  that  the  plans 
computed  by  what  are  called  optimization  systems  are  not 
consistent  with  the  expectation  of  the  planner  and  that  sev¬ 
eral  trial-and-error  adjustments  of  the  system  parameters 
might  be  required  to  achieve  a  clinically  acceptable  plan. 
Given  a  patient,  the  obtained  plan  can  vary  widely  from  one 
planer  to  the  next,  even  within  a  department.  In  the  following 
we  describe  a  more  adaptable  and  “intelligent”  statistical 
inverse  planning  formalism  based  on  the  concept  of  a  pref¬ 
erence  function  to  better  deal  with  the  dilemma. 

B.  Preference  function 

In  a  recent  paper,  Xing  et  a/.31  introduced  the  concept  of 
preference  function  to  weaken  the  rigid  dose  prescription 
commonly  seen  in  the  existing  inverse  planning  algorithms. 
Its  role  is  to  allow  a  dose  distribution  to  be  considered  in¬ 
stead  of  just  a  single  value,  and  to  quantify  the  degree  of  our 
willingness  to  accept  a  prescription  dose  Dp  in  that  range. 
The  preference  function  can  be  constructed  heuristically 
from  clinical  considerations.31  The  defined  preference  func¬ 
tion  states  that  the  most  favorable  prescription  dose  for  a 
voxel  n  is  Dp(n)  and  that  a  different  prescription  dose  is  also 
acceptable,  but  with  a  smaller  preference  level.  For  illustra¬ 
tion,  in  Fig.  1  we  show  a  sketch  of  the  preference  functions 


for  a  target  and  sensitive  structure.  The  most  desirable  dose 
for  a  sensitive  structure  should  generally  be  set  to  zero.  The 
conventional  prescription  scheme  represents  a  special  case  of 
the  general  approach  proposed  here  with  the  step  function 
form  of  the  preference  function.  That  is, 


( 1,  if  Dp  =  D°p, 

P"(D')“|0,  if  Dp±D°p. 


(3) 


To  give  another  example,  we  write  down  the  Gaussian 
preference  function  for  a  voxel  n 

P^Dp)  =  P0„exp{-y„[Dp(n)-D°p(n)]2},  (4) 


where  P0n  is  a  normalization  constant  and  y„  represents  the 
Gaussian  parameter.  For  a  system  comprising  N  voxels,  the 
total  preference  is  given  by  a  product  of  the  preference  func¬ 
tions  of  all  voxels: 


p=u  p*{dp) 

n 


=  11  Po^w{-yn[Dp{n)-D°p{n)f}.  (5) 

n 

When  a  maximum  likelihood  estimator  is  used,  it  has 
been  demonstrated  that  the  maximization  of  the  logarithmic 
function  of  P  or  minimization  of  lnjl/P),  is  equivalent  to  the 
minimization  of  the  conventional  quadratic  objective 
function.31,34  In  this  case,  the  Gaussian  parameter  y„  in  Eq. 
(5),  which  commands  the  “spread”  of  the  Gaussian  around 
Dp,  is  equivalent  to  the  importance  factor  that  controls  the 
relative  importance  of  the  structure  and  parametrizes  the 
clinical  trade-off  strategy. 

C.  Probability  density-based  dose  prescription  and 
inverse  planning 

The  objective  function  defined  in  Eq.  (2)  uses  a  rigid 
dose,  Dp .  Since  in  most  instances  an  ideal  dose  prescription 
is  not  physically  attainable,  we  resort  to  an  expansion  of  the 
prescription  dose,  over  a  certain  interval.  That  is,  we  allow 
the  prescription  dose  to  take  a  “probabilistic”  distribution 
around  the  most  desired  dose  as  specified  by  the  preference 
function.  For  computational  purpose,  we  divide  the  permis- 
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sible  prescription  dose  into  a  number  of  discretized  values, 
{D'p},  where  i  is  the  index  of  a  possible  prescription  dose  and 
i  =  0  represents  the  most  desirable  dose.  The  preference  dis¬ 
tribution  prescription  is  usually  normalized  to  unity. 

In  order  to  utilize  the  probability  information  character¬ 
ized  by  the  preference  function,  we  formulate  the  conven¬ 
tional  dose  optimization  into  a  statistical  analysis  problem. 
To  proceed,  let  us  take  the  quadratic  objective  function  as  an 
example.  We  rewrite  the  traditional  quadratic  objective  func¬ 
tion  (2)  into 

/(DC)=/0][I  exp {-ra[Dc{n)-Dp{n)']2}.  (6) 

n 

where  /0  is  a  normalization  constant.  For  a  given  prescribed 
dose  distribution,  Eq.  (6)  measures  the  goodness  of  a  calcu¬ 
lated  dose  distribution  using  an  exponential  scale,  as  com¬ 
pared  with  Eq.  (2).  Equation  (6)  can  be  interpreted  as  a  con¬ 
ditional  probability  and  formally  rewritten  as 

f{Dc\Dp)  =/0II  exp {-ra[Dc{n)-Dp{n)~\2}.  (7) 

n 

When  the  prescription  dose  is  no  longer  a  rigid  dose,  it  is 
conceivable  that  there  are  a  number  of  optimum  solutions, 
each  corresponding  to  a  sample  of  prescription  doses.  Math¬ 
ematically,  we  now  have  two  “probability”  distribution  func¬ 
tions.  One  is  the  preference  function  that  characterizes  our  a 
priori  preference  over  different  prescription  doses  P(Dp), 
and  the  other  is  Eq.  (6)  that  ranks  a  calculated  dose  for  a 
given  prescribed  dose,  Dp .  Our  task  is  to  find  the  solution 
that  is  statistically  optimal  with  consideration  of  the  variable 
prescription.  For  this  purpose,  we  introduce  the  “joint  prob¬ 
ability”  of  the  two  “probability”  distributions  defined  by 
Eqs.  (5)  and  (7).  The  function  at  a  voxel  n  can  be  written  as 

P„(DC)^  UD^D^P^).  (8) 

/ 

The  total  preference  function  of  the  system  is  given  by 

^=n  pn(Dc).  (9) 


D.  Optimization  strategy 

Having  the  rigid  prescription  Dp  in  (2)  replaced  by  a 
range  of  prescribed  doses,  {D'p},  the  total  preference  func¬ 
tion  is  now  given  by  Eqs.  (8)  and  (9).  For  convenience,  we 
define  objective  function  F=ln(l/P)  and  derive  the  optimal 
solution  by  minimizing  the  F,  which  is  equivalent  to  maxi¬ 
mize  the  preference  function  (9).  The  objective  function  now 
reads  as 

F=\n(\/P) 

=  -lnIl^=-S  'n2  f n(D c\D'p) ■  P „(D'p) .  (10) 

n  n  i 


Fig.  2.  A  flow  chart  of  the  optimization  process  with  the  inclusion  of  pre¬ 
designed  preference  function  information. 

Note  that  the  conventional  quadratic  objective  function  is  a 
special  case  of  the  above  general  objective  function  when  the 
prescription  takes  a  rigid  value  for  each  structure,  as  de¬ 
scribed  by  Eq.  (3). 

The  optimization  process  is  schematically  shown  in  a  flow 
chart  (Fig.  2).  The  beam  profile  is  determined  by  minimizing 
the  above  objective  function  using  a  conjugate  gradient  op¬ 
timization  algorithm.  The  details  of  the  algorithm  have  been 
discussed  in  a  previous  paper.33  Briefly,  the  calculation  con¬ 
sists  of  three  major  steps:  (i)  assume  an  initial  intensity  pro¬ 
file  for  each  incident  beam;  (ii)  compute  the  “joint  probabil¬ 
ity”  given  by  Eqs.  (8)  and  (9).  For  this  purpose,  we  need  to 
sample  all  combinations  of  the  prescription  doses  of  different 
structures  and  compute  the  function  given  in  Eqs.  (5)  and  (7) 
for  each  of  these  combinations;  and  (iii)  optimization  of  the 
multidimensional  “joint  probability”  function.  The  second 
step  is  fairly  computationally  intensive  because  we  must 
compute  the  two  functions  for  every  sampling  of  the  pre¬ 
scription  doses.  In  our  calculation,  we  typically  assign  four 
to  seven  discrete  possible  prescription  doses  for  each  struc¬ 
ture.  A  finer  discretization  of  the  prescription  dose  did  not 
seem  lead  to  further  improvement  but  would  greatly  increase 
the  computation  time.  All  calculations  presented  here  are 
performed  on  a  Personal  Computer  (PC)  with  an  Intel  Pen¬ 
tium®  III  1  GHz  CPU  (Intel  Corporation,  Sunnyvale,  CA). 
The  computation  time  needed  to  obtain  an  optimal  solution 
for  a  given  set  of  system  parameters  (including  beam  con¬ 
figuration,  preference  function,  importance  factors)  is  typi¬ 
cally  less  than  ten  minutes. 

III.  RESULTS  AND  DISCUSSION 

A.  A  synthetic  phantom  case  with  a  C-shaped  tumor 

To  systematically  study  the  performance  of  the  statistical 
analysis-based  inverse  planning  algorithm,  we  applied  the 
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Fig.  3.  (a)  A  sketch  of  a  phantom  case  with  a  C-shaped 
tumor.  The  dose  prescription  is  set  100  (arbitrary  units)  to 
the  PTV  and  0  to  the  circular  OAR  and  normal  tissue,  (b) 
Dose  distribution  obtained  using  the  “probabilistic”  pre¬ 
scription  shown  in  Fig.  4(a). 


technique  to  a  C-shaped  tumor  case  [Fig.  3(a)]  with  a  variety 
of  preference  functions  and  compared  the  results  with  that 
obtained  using  the  conventional  approach  with  a  fixed  dose 
prescription.  Nine  equally  spaced  6  MV  beams  beginning  at 
0°  (IEC)  were  used  in  this  study.  The  prescription  doses  to 
the  PTV  and  OAR  in  the  conventional  IMRT  plans  were  100 
and  0  (the  dose  is  in  an  arbitrary  unit),  respectively. 

We  first  assigned  three  sets  of  symmetrical  Gaussian  dis¬ 
tributions  to  the  target  while  keeping  the  prescription  to  the 
sensitive  structure  at  zero  (Fig.  4).  The  Gaussian  preference 
functions  were  represented  by  three  sets  of  preference  levels 
at  seven  discrete  values  (80,  87,  94,  100,  106,  113,  and  120). 
The  center  of  the  Gaussian  functions  was  set  at  100.  The 
preference  levels  for  the  seven  doses  are  shown  in  Fig.  4  for 
each  of  the  three  situations  studied  here.  The  transverse  dose 
distribution  obtained  using  the  statistical  inverse  planning 
formalism  for  the  case  shown  in  Fig.  4(a)  is  plotted  in  Fig. 
3(b).  As  expected,  target  inhomogeneity  increases  as  we 


loosen  the  constraint  of  the  rigid  dose  prescription.  This  can 
be  better  demonstrated  by  using  the  differential  DVH  for 
each  situation.  As  seen  from  the  differential  DVH  plots  (the 
right  column  of  Fig.  4),  the  width  of  the  differential  function 
gradually  increases,  from  26.72,  28.59-30.39,  as  we  gradu¬ 
ally  increase  the  acceptance  levels  for  the  doses  different 
from  the  most  desirable  dose  (100).  This  series  of  calcula¬ 
tions  provides  us  with  preliminary  evidence  that  the  final 
dose  distribution  can  be  steered  by  varying  the  preference 
function. 

Next,  we  constructed  six  sets  of  asymmetric  preference 
functions  for  the  target  (Fig.  5  and  Fig.  6).  When  higher 
preference  levels  were  assigned  to  the  doses  higher  than  100, 
we  found  that  the  target  DVH  is  shifted  to  the  high  dose 
region.  Interestingly,  even  when  an  extremely  low  preference 
(for  instance,  1%)  was  assigned  to  the  doses  less  than  100 
[Fig.  5(b)],  a  noticeable  underdosing  relative  to  the  conven¬ 
tional  result  was  resulted.  A  similar  phenomenon  can  also  be 
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Dose  (atitrary  unit)  Dose(attiayirt) 

Dose  (arbitrary  unit) 


B 


Dose  (arbitrary  unit) 


Dose(abitiayuit)  Dose(abrtrayuit) 


Dose  (arbitrary  unit) 


Fig.  4.  DVHs  of  the  PTV,  OAR,  and  normal  tissue  (NT)  obtained  using  the  conventional  rigid  dose  prescription  (dotted  line)  and  the  “probabilistic” 
prescription  (solid  line).  The  Gaussian  preference  functions  with  different  variances  are  shown  in  the  middle  panel.  The  right  panel  shows  the  differential 
DVHs  for  the  target. 
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... —  uniform 
- Probabilistic 
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Dose  (arbitrary  irrt) 
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Dose  (arbitrary  unit) 


as- 


SO  87  94  100  108  113  130 


Dose  (arbitrary  irit) 


Fig.  5.  DVHs  of  the  PTV,  OAR,  and  normal  tissue  (NT)  obtained  with  the  conventional  rigid  dose  prescription  (dotted  line)  and  with  the  “probabilistic” 
prescription  (solid  line).  The  bar  charts  on  right  show  the  asymmetrical  preference  functions. 
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0  20  4b  60  80  100  120  DosefaifoHi  <iiy  unit) 


Dose  (arbitrary  unit) 
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Dose  (arbitrary  unit) 


80  87  94  100  «6  113  120 

Dose  (arbitrary  unit) 


Dose  (arbitrary  unit)  Dose(atttrayint) 

Fig.  6.  DVHs  of  the  PTV,  OAR,  and  normal  tissue  (NT)  obtained  using  the  conventional  rigid  dose  prescription  (dotted  line)  and  the  new  statistical  inverse 
planning  method  for  a  variety  of  preference  functions  shown  on  the  right  panel  (solid  line). 


seen  from  the  result  shown  in  Fig.  5(c),  where  only  0.5%, 
0.7%,  and  1%  of  preference  levels  were  assigned  to  the  dose 
values  of  80,  87,  and  94.  This  observation  seems  to  indicate 
that  the  influence  of  the  assigned  preference  level  at  a  low 


dose  plays  an  important  role.  In  Fig.  6,  we  set  the  preference 
levels  for  the  doses  less  than  100  to  be  0  and  only  assign 
nonzero  preference  levels  for  the  doses  higher  than  100.  It  is 
seen  that  in  all  these  situations  the  minimum  target  dose  is 
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OAR  Preference 


(A)  <B) 


0  5  10  15  20  25  90 

Dose  (arbitrary  unit) 


(D) 


Fig.  7.  DVHs  of  the  OAR  and  PTV  when  the  prescription  dose  to  the  OAR  is  modeled  by  (a)  a  uniform  distribution,  (b)  a  bell-shaped  function,  (c)  an 
exponential  decay  function,  and  (d)  a  rigid  value. 


higher  than  that  of  the  conventional  plan.  As  a  result  of  our 
preference  over  higher  doses,  the  fractional  volume  at  any 
dose  less  than  100  is  improved  in  comparison  to  that  of  the 
conventional  IMRT  plan.  In  Fig.  6(c),  we  further  exemplify 
the  statistical  analysis  based  inverse  planning  method  by 
simplifying  our  preference  to  two  doses  (100  and  120),  each 
with  50%  preference  levels.  In  this  situation,  in  addition  to 
that  the  doses  in  the  target  are  shifted  toward  higher  values, 
the  target  DVH  exhibits  a  stepwise  behavior:  a  plateau  ap¬ 
pears  at  around  1 10,  which  is  in  the  middle  of  the  two  pre¬ 
scribed  doses. 

It  is  interesting  to  point  out  that  the  OAR  sparing  is  im¬ 
proved  as  compared  with  the  conventional  IMRT  plan  in 
most  cases  studied  in  Figs.  5  and  6,  even  when  the  target 
dose  is  escalated.  That  is,  the  DVH  of  the  OAR  is  not  always 
shifted  toward  higher  doses,  as  would  occur  if  a  higher  dose 
is  prescribed  in  a  conventional  inverse  planning  system.  In¬ 
stead,  the  dose  to  the  OAR  remained  unchanged  or  even 
lowered  in  some  cases.  A  reasonable  explanation  for  the  ob¬ 
served  phenomenon  is  that,  when  a  rigid  dose  prescription  is 
replaced  by  a  range  of  doses,  the  system  is  given  more  free¬ 
dom  for  self-adjustment.  As  a  benefit,  a  solution  with  a 
higher  integral  target  dose  and  reduced  OAR  dose  can  be 
obtained  from  the  expanded  solution  space. 

We  have  also  studied  the  behavior  of  the  system  when  a 
range  of  doses  is  prescribed  to  the  OAR.  In  this  investiga¬ 
tion,  we  kept  the  target  prescription  to  100  and  allowed  the 


OAR  dose  to  take  seven  values:  0,  5,  10,  15,  20,  25,  and  30 
with  the  acceptance  levels  sampled  from  three  different  types 
of  prescription  distribution:  uniform  [Fig.  7(a)],  bell-shaped 
[Fig.  7(b)],  and  exponential  [Fig.  7(c)]  functions.  Figure  7(d) 
represents  the  conventional  case  with  zero  prescription  to  the 
OAR.  The  corresponding  OAR  and  PTV  DVHs  are  plotted 
in  the  left  panel  of  Fig.  7.  When  the  preference  was  uni¬ 
formly  sampled  in  the  dose  interval  from  0  to  30,  the  result¬ 
ant  dose  to  the  OAR  was  found  to  be  the  highest,  as  indi¬ 
cated  by  curve  A  in  Fig.  7.  The  best  target  dose  coverage  was 
achieved  in  this  situation.  If  the  preference  to  a  high  dose 
was  reduced,  the  DVH  was  gradually  shifted  to  the  low  dose 
direction  (curves  B  and  C).  It  is  not  surprising  that  the  best 
OAR  sparing  was  achieved  in  the  conventional  case  where  a 
zero  dose  was  prescribed  to  the  OAR.  The  target  dose  homo¬ 
geneity  was  slightly  improved  in  all  cases  when  a  probabi¬ 
listic  prescription  was  given  to  the  OAR.  Similar  to  that  de¬ 
scribed  in  the  last  paragraph,  the  results  clearly  demonstrate 
that  the  “probabilistic”  prescription  allows  us  to  control  the 
OAR  dose  distribution  and  indicate  the  usefulness  of  the 
statistic  analysis  approach. 

B.  The  prostate  case 

The  new  inverse  planning  algorithm  was  also  applied  to 
study  a  six  filed  IMRT  prostate  treatment  [Fig.  8(a)].  Four 
plans  with  different  types  of  preference  functions  were  gen- 
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Fig.  8.  A  transverse  slice  showing  the  anatomical  struc¬ 
tures  delineated  for  the  prostate  tumor  (a)  and  (b)  the 
dose  distribution  obtained  using  the  “probabilistic” 
prescription  shown  in  Fig.  9(a). 


erated.  In  addition,  a  plan  with  rigid  prescription  (74  Gy  on 
the  target,  60  Gy  on  the  bladder,  and  40  Gy  on  the  rectum)  is 
also  generated.  The  DVHs  for  this  plan  is  plotted  as  dotted 
lines  in  Fig.  9  and  is  used  as  a  reference  for  comparison.  In 
all  treatment  plans,  six  beams  were  placed  at  the  following 
angular  positions:  0°,  55°,  135°,  180°,  225°,  and  305°.  The 
size  of  the  pencil  beam  defined  at  the  isocenter  was  0.5  cm. 

The  DVFI  and  preference  functions  for  four  different 
plans  are  schematically  shown  in  Fig.  9.  In  the  study  shown 
in  Figs.  9(a)-9(b),  we  kept  the  preference  function  of  the 
sensitive  structures  unchanged  and  only  varied  the  form  of 
the  preference  function  of  the  target.  In  Fig.  9(a),  we  as¬ 
sumed  that  target  could  take  seven  discrete  values  (74,  76, 
78,  80,  82,  84,  and  86  Gy)  sampled  from  an  exponential 
distribution.  Compared  with  the  dotted  lines,  the  target  DVH 
was  shifted  toward  the  high  dose  direction.  The  dose  distri¬ 
bution  corresponding  to  the  preference  function  is  shown  in 
Fig.  8(b).  The  target  DVH  was  shifted  even  further  toward 


the  high  dose  region  [Fig.  9(b)]  when  a  bell-shaped  prefer¬ 
ence  function  was  used  with  more  emphasis  on  the  target 
receiving  doses  at  74,  76,  and  78  Gy.  In  both  cases,  doses  to 
the  rectum  and  bladder  did  not  change  significantly. 

In  Fig.  9(c)  we  show  the  DVHs  when  the  preference  func¬ 
tion  to  the  rectum  deviates  from  the  uniform  distribution.  As 
a  result,  the  rectum  dose  was  significantly  lowered  in  all  dose 
levels  and  the  maximum  dose  was  reduced  from  66  to  57  Gy. 
Because  of  the  proximity  of  the  rectum  to  the  prostate  target, 
the  maximum  rectum  dose  was  not  restricted  to  30  Gy,  as 
specified  in  the  preference  function.  We  emphasize  that  the 
improvement  in  rectum  and  bladder  sparing  was  achieved  at 
cost  of  higher  dose  inhomogeneity  in  the  prostate  target.  This 
reminds  us  that,  in  dose  optimization,  there  is  a  dosimetric 
compromise.  That  is,  the  improvement  in  the  dose  to  a  struc¬ 
ture  is  often  accompanied  by  dosimetrically  adverse  effect(s) 
at  other  points  in  the  same  or  different  structures.  The  impor¬ 
tant  point  that  one  should  note  is  that  from  the  clinical  point 
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Fig.  9.  A  comparison  of  the  DVHs  obtained  using  the  “probabilistic”  prescription  (solid  line)  and  conventional  rigid  dose  prescription  (dotted  line).  The 
prostate  and  rectum  prescriptions  are  represented  on  the  right  bar  charts. 
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of  view,  some  dose  distributions  are  more  acceptable  than 
others  and  our  goal  is  to  find  the  solution  that  improves  the 
plan  to  the  largest  possible  extent,  but  with  a  clinically  insig¬ 
nificant  or  acceptable  sacrifice  in  OAR  sparing.  In  order  to 
achieve  this,  it  is  necessary  to  have  a  reasonable  amount  of 
controllability  over  the  final  dose  distribution.  In  this  sense 
we  believe  that  the  proposed  formalism  is  valuable. 

In  addition,  we  varied  the  preference  functions  for  both 
target  and  rectum  [Fig.  9(d)].  The  preference  function  for  the 
rectum  was  the  same  as  that  in  Fig.  9(c).  Compared  to  the 
results  shown  in  Fig.  9(c),  we  found  that  the  target  dose 
inhomogeneity  was  slightly  improved. 


IV.  CONCLUSIONS 

The  formalism  we  derived  here  provides  a  general  starting 
point  for  the  study  of  a  system  with  a  probability  density- 
based  dose  prescription.  The  inclusion  of  the  partial  informa¬ 
tion  into  the  plan  selection  process  represents  a  significant 
change  from  the  conventional  approaches.  The  proposed 
technique  can  be  categorized  into  the  general  Bayesian 
decision-making  theory,31  which  is  a  useful  tool  to  deal  with 
a  system  with  “statistical”  inference.  In  image  analysis  and 
many  other  fields  of  science  and  engineering,  it  has  proven 
extremely  useful  to  include  the  prior  knowledge  of  the  sys¬ 
tem  variables  into  the  estimation  process.31  The  preference 
function  proposed  for  radiotherapy  optimization  here  serves 
as  a  priori  probability  density  function  in  standard  Bayesian 
statistics.  The  role  of  the  preference  function  is  to  indicate 
our  “bias”  on  the  values  of  the  system  variables.  By  utilizing 
the  partial  information  of  the  system  variables,  one  can  more 
effectively  search  the  solution  space  and  eliminate  some  un¬ 
necessary  uncertainties  in  the  optimization  process. 

In  conclusion,  we  have  developed  a  statistical  analysis- 
based  inverse  planning  algorithm  to  include  preference  and 
expert  knowledge  into  the  dose  optimization  process.  Instead 
of  a  rigid  dose  prescription,  the  new  approach  allows  us  to 
prescribe  a  range  of  doses  with  predesigned  preference  lev¬ 
els.  The  technique  represents  a  novel  application  of  the  gen¬ 
eral  Bayesian  decision-making  theory31  for  dealing  with  sta¬ 
tistical  inference  and  is  valuable  for  deriving  a  statistically 
optimal  solution  in  the  presence  of  uncertainties  in  system 
parameters.  The  method  was  demonstrated  for  a  system  with 
modulating  prescriptions  but  can  be  easily  extended  to  solve 
many  other  related  problems  (e.g.,  in  biologically  based  dose 
optimization,  one  can  incorporate  the  uncertainties  of  various 
radiobiology  parameters  into  the  inverse  planning  process 
using  the  frameset  developed  in  this  work35).  The  ill  condi¬ 
tioning  of  the  problem  was  improved  because  of  the  use  of  a 
less  restrictive  prescription  and,  as  a  result,  new  solutions 
that  are  otherwise  inaccessible  can  be  obtained  naturally.  It  is 
demonstrated  that  the  obtained  solutions  using  the  new  ap¬ 
proach  strongly  correlate  with  the  preference  function,  sug¬ 
gesting  that  the  planning  process  is  controllable  and  predict¬ 
able  by  the  proposed  method. 
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